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Ultra-Low  Threshold  Vertical-Cavity  Surface-Emitting  Lasers  for 

USAF  Applications 

Thomas  R.  Nelson,  Jr.,  PhD. 

AFRL/SNDD,  2241  Avionics  Circle, 

Wright  Patterson  AFB,  OH  4 54 HP 

(Dated:  January  11,  2005.) 

Abstract 

This  report  summarizes  research  on  the  development  of  ultra-low  threshold  vertical-cavity 
surface-emitting  lasers.  This  in-house  lab  task,  supported  by  AFOSR  under  LRIR  2305DW01, 
was  initiated  by  Dr.  John  Loehr,  with  AFOSR  Program  Manager  Dr.  Alan  Craig.  In  2001,  Dr. 
Thomas  Nelson  became  the  in-house  project  lead  on  this  effort,  with  AFOSR  management  trans¬ 
ferred  first  to  Dr.  Kent  Miller,  and  most  recently  to  Dr.  Gernot  Pomrenke.  This  task  was  initiated 
in  March  1996,  and  this  report  covers  roughly  the  time  period  of  March  1996  -  Dec  2002. 
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INTRODUCTION:  GOAL  OF  THIS  RESEARCH 


The  goal  of  this  Air  Force  Office  of  Scientific  Research  (AFOSR)  sponsored  in-house  labo¬ 
ratory  task  has  been  the  development  of  ultra-low  threshold  vertical-cavity  surface  emitting 
lasers  (VCSELs)  for  USAF  applications.  Interest  in  VCSELs  by  the  Air  Force  Research 
Laboratory  (AFRL)  Sensors  Directorate,  specifically  from  the  Aerospace  Components  and 
Subsystems  Concept  Division,  Electron  Devices  Branch  (AFRL/SNDD)  was  in  several  are¬ 
nas: 

*  Low  cost,  short  haul  optical  local  area  networks  (LANS),  with  the  vision  of  a  ” fly-by- 
fiber”  approach  to  future  aircraft  systems. 

*  Phased  arrays  of  VCSELs  for  potential  laser  radar  (LADAR)  applications. 

*  High  speed  optical  components  for  implementation  in  radio-frequency  to  photonic 
links  (RF /Photonic  links)  for  cost  effective  information  transmission  of  conventional 
RF  radar  signals. 

*  The  development  of  photonic  elements  involved  in  optical  true-time-delay  systems, 
again  for  improving  conventional  RF  radar  systems  by  the  use  of  optoelectronic  sub¬ 
systems. 

In  all  of  these  efforts,  several  underlying  goals  were  desirable.  These  include  low-cost 
mass  production  methods,  high-speed  device  operation,  ease  of  integration  onto  existing 
or  future  platforms,  and  a  robust  device  capable  of  withstanding  the  unique  environments 
experienced  by  air  and  space  craft. 

At  the  time  of  the  initiation  of  this  in-house  laboratory  task,  semiconductor  optoelec¬ 
tronic  devices,  including  laser  diodes,  were  advancing  from  fundamental  research  concepts 
to  advanced  engineering  projects  for  inclusion  in  photonic  components  and  subsystems.  In¬ 
deed,  a  few  companies  were  emerging  that  for  the  first  time  were  making  vertical-cavity 
surface-emitting  lasers  (VCSELs)  commercially  available.  Early  emphasis  was  on  high- 
bandwidth  applications  and  developing  arrays  of  such  devices  for  transceiver  modules  in 
telecommunications  systems.  Such  products  were  based  on  high-yield  processes,  not  nec¬ 
essarily  state-of-the-art  optoelectronics  fabrication  methods  designed  at  high-performance 
devices.  Furthermore,  as  their  usage  was  intended  for  terrestrial  systems,  overall  power 
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consumption,  heat  dissipation,  etc  were  not  effects  of  enormous  concern.  On  an  aerospace 
platform,  however,  total  power  budget  is  a  primary  driving  force  in  total  system  design,  and 
the  use  of  high  efficiency  components  is  mandatory.  These  concerns,  among  others,  gave  the 
uitra-iow  threshold  focus  to  the  effort  summarized  in  this  document. 

In  ail  the  work  detailed  below,  unless  otherwise  noted,  the  research  was  performed  entirely 
“in-house.”  To  clarify,  the  effort  started  with  theoretical  device  design  based  upon  the 
work  of  AFRL/SNDD  scientists,  working  in  conjunction  with  the  Air  Force  Institute  of 
Technology  (co- located  at  Wright-Patterson  Air  Force  Base).  Device  material  growth  was 
accomplished  by  means  of  molecular  beam  epitaxy  (MBE)  using  a  Varian  GENII  MBE 
system  owned  by  the  Materials  Directorate  (AFRL/MLPSM).  This  growth  system  is  located 
in  the  cleanroom  of  the  Sensors  Directorate,  is  operated  by  AFRL/MLPSM  scientists,  and 
many  upgrades  and  modifications  were  made  to  this  tool  to  improve  its  device  growth 
capability  as  a  direct  result  of  the  stringent  growth  requirements  of  this  project.  These  will 
be  detailed  in  later  sections.  Similarly,  all  device  fabrication  and  characterization  were  also 
accomplished  at  AFRL/SNDD.  One  significant  accomplishment  as  a  result  of  this  effort  was 
the  development  of  a  steam  oxidation  furnace  with  in-situ  oxidation  monitoring,  and  its  role 
in  device  fabrication  will  be  explained  in  later  sections. 

BACKGROUND  AND  THEORY 

This  project  officially  commenced  in  March  1996,  under  the  guidance  of  in-house  project 
leader  Dr.  John  Loehr,  and  AFOSR  program  manager  Dr.  Alan  Craig.  The  first  goal 
was  to  determine  exactly  what  type  of  VCSEL  structure  to  pursue.  An  overview  VCSEL 
taxonomy  is  presented  in  Fig.  1  showing  the  most  common  three  varieties  of  such  devices, 
namely  the  etched  post  VCSEL  (Fig.  1A),  the  proton-implanted  VCSEL  (Fig.  IB),  and 
the  intra-cavity  contacted  VCSEL  with  native-oxide  aperture  layers  (Fig.  1C).  From  these 
examples,  several  design  issues  are  apparent:  VCSELs  can  be  completely  monolithic  in 
nature,  with  both  top  and  bottom  distributed  Bragg  reflector  (DBR)  mirrors  formed  from 
epitaxial  growth  (as  in  Fig.  1A  and  Fig.  IB),  or  partially  monolithic  using  a  combination 
of  epitaxial  (bottom  mirror,  Fig.  1C)  and  post-growth  deposition  (top  mirror,  Fig.  1C); 
VCSELs  can  have  their  current  contacts  formed  outside  the  cavity  by  suitably  doping  the 
DBR  mirrors  (Fig.  1A  and  Fig.  IB),  or  via  intra-cavity  contacting  (Fig.  1C);  and  there  are 
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a  variety  of  means  for  confining  the  optical  field  as  well  as  the  injected  current  via  etched 
posts,  proton  implantation  followed  by  annealing,  native  oxide  apertures,  or  selective  etching 
(not  shown).  In  each  case,  significant  tradeoffs  exist.  For  instance,  proton  implanted  devices 
without  mesa  isolation  etching  received  the  earliest  commercial  interest  due  to  their  relative 
ease  of  fabrication.  Production  of  these  devices  requires  extensive  initial  calibrations  of  the 
proton  implant  profile,  however,  not  to  mention  the  availability  of  an  implantation  system. 
Etched  post  structures  provide  better  optical  confinement,  but  the  deep  etching  through 
the  DBR  mirrors  required  in  these  structures  warrants  exploring  the  relative  benefits  of  wet 
versus  dry  etching,  in  addition  to  the  potential  need  to  passivate  the  sidewalls  roughened 
by  such  etching  to  avoid  surface  recombination  losses.  Finally,  the  native-oxide  aperture 
devices  require  not  only  the  etched  post  techniques  above,  but  also  an  apparatus  (steam 
bubbler  type  system)  running  at  sufficiently  high  temperature  (T  >  375 °C)  to  convert  high 
mole  fraction  Al(.i')Ga(  1  —  a;) As  (x  >  0.95)  material  into  native  oxides  of  aluminum. 

In  order  to  determine  which  structure  to  study,  a  program  was  initiated  to  develop  a  suite 
of  modelling  tools  to  allow  us  to  compare  the  relative  benefits  of  each  structure.  Much  of 
the  bandstructure  code  useful  in  determining  laser  gain  as  a  function  of  injected  carriers  was 
already  in  place  from  previous  research  [1-3],  and  could  readily  be  incorporated  into  this 
suite.  More  challenging,  however,  was  the  ability  to  model  the  optical  “cold  cavity”  fields 
in  structures  such  as  the  oxide- apertured  etched-post  designs.  The  goal  in  this  design  effort 
is  the  reduction  (ideally,  elimination)  of  threshold  current,  that  is,  carriers  (electrons  and 
holes)  that  recombine  but  do  not  contribute  to  lasing.  In  large  area  lasers  (edge-emitting 
or  VCSELs),  the  threshold  current  is  proportional  to  the  volume  of  the  active  region.  In 
most  designs,  this  active  region  consists  of  one  or  more  quantum  wells  (QWs)  separated  by 
appropriate  thickness  barriers,  and  multiple  sets  of  such  wells  and  barriers  are  then  grown 
at  appropriate  positions  in  the  lasing  cavity  so  as  to  have  maximum  spatial  overlap  with 
the  intracavity  lasing  field  mode.  The  use  of  both  etched-post  and  oxide- apertured  struc¬ 
tures  in  early  efforts  demonstrated  improvements  in  threshold  current  reduction,  primarily 
by  restricting  and  refining  the  transverse  gain  profile  to  overlap  better  with  the  optical  field 
mode.  However,  worth  noting  is  the  fact  that  such  lateral  patterning  also  affects  the  trans¬ 
verse  optical  confinement  factor.  When  the  posts  and/or  apertures  shrink  to  dimensions  on 
the  order  of  the  lasing  emission  wavelength,  three  optical  effects  begin  to  emerge:  (1)  The 
spatial  profile  of  each  cavity  mode  changes  from  a  plane  wave  to  a  3-D  waveguide  type  of 
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FIG.  1:  Schematic  cross  sections  of  three  common  types  of  VCSELs.  (A)  An  etched  post  structure, 
extra-cavity  contacted,  with  oxide  apertures  (B)  A  proton  implanted,  extra-cavity  contacted  device. 
(C)  An  intra-cavity  contacted  device  with  oxide  apertures,  showing  the  possibility  of  top  and 
bottom  DBR  mirrors  formed  from  different  methods  (epitaxially  grown  bottom  mirror,  deposited 
top  mirror). 

mode;  such  a  change  then  alters  the  modal  reflectivity  and  the  optical/active  region  overlap; 
(2)  The  boundaries  of  the  oxide  aperture  or  etched  post  begin  to  act  as  sources  of  diffraction 
losses,  shedding  optical  radiation  laterally  from  the  device;  (3)  finally,  the  cavity  lasing  mode 
begins  to  blueshift  in  wavelength,  with  the  spacing  between  modes  increasing.  All  of  these 
effects,  then,  play  an  important  role  in  the  effective  design  of  low  threshold  lasers. 

Previous  efforts  to  simulate  VCSEL  cavity  structures  incorporating  oxide  apertures  have 
varied  significantly  in  their  level  of  complexity  and  detail,  in  addition  to  their  predicted 
results.  Some  groups  suggested  that  thin  oxides  placed  at  field  nodes  provide  the  lowest 
threshold,  while  others  argued  that  thick  oxides  at  antinodes  are  the  correct  approach.  The 
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disparity  in  answers  is  largely  due  to  an  incomplete  understanding  of  the  optical  physics 
related  to  the  aperture,  or,  more  precisely,  the  aperture/cavity  optical  system.  In  an  attempt 
to  understand  the  physics,  several  models  have  been  created,  all  with  different  approaches 
to  the  problem.  However  none  of  these  models  provide  a  comprehensive  description  of  the 
optical  fields,  and  no  clear,  consistent  design  guidance  existed.  In  order  to  overcome  these 
limitations,  we  embarked  on  efforts  to  develop  methods  for  modelling  lasing  modes  of  various 
VCSEL  structures. 

Two  sets  of  simulation  codes  were  developed  to  aid  in  ultra-low  threshold  VCSEL  design. 
The  first  of  these,  the  Weighted  Index  Method  with  Parasitics  (WIMP),  is  a  semi-analytic 
hybrid  vector-scalar  code.  As  such,  it  was  designed  to  be  implemented  on  relatively  standard 
computing  platforms.  Indeed,  the  source  code  that  our  group  uses  was  written  in  Visual 
Fortran  for  the  Microsoft  Windows  (XP,  ME,  or  9x)  platforms,  and  it  doesn’t  require  unrea¬ 
sonably  heavy  duty  computing  power  (large  RAM  or  CPU  requirements,  for  example).  The 
only  inflexible  portions  of  the  code  at  this  point  in  time  are  the  usage  of  bundled  numerical 
library  routines  from  the  International  Mathematical  and  Statistical  Libraries  (IMSL)  suite 
of  software  codes.  These  library  routines  perform  such  functions  as  complex  root  finding, 
calculation  of  Bessel  functions,  and  numerical  integration  and  differentiation.  Typical  run 
times  for  problems  seeking  to  find  the  first  few  roots  of  an  oxide  aperture  VCSEL  of  typical 
dimensions  (etched  post  of  ~  20/un  diameter,  oxide  aperture  diameter  of  ~  5 //in)  usually 
run  less  than  10  minutes  on  a  Pentium  4  1-GHz  CPU  platform.  Drawbacks  and  limitations 
of  this  method  will  be  explained  in  the  later  section  detailing  simulation  results. 

The  second  method  used  for  simulation  of  etched-post  and  oxide  aperture  VCSELs  is 
the  Vector  Finite  Element  Method  (VFEM).  Unlike  WIMP,  this  is  a  numerically  intensive 
method  due  to  the  meshing  of  potentially  hundreds  of  layers  forming  an  oxide  aperture  VC¬ 
SEL,  and  also  due  to  the  fact  that  the  basis  elements  are  vector  and  not  scalar  in  nature. 
Indeed,  as  will  be  shown  later,  a  set  of  14  basis  vector  functions  (6  node  based,  6  edge  based, 
and  2  face  based)  were  chosen  for  this  model.  Physical  RAM  computing  requirements  dic¬ 
tated  that  these  simulations  be  performed  at  the  Major  Shared  Resources  Center  (MSRC),  a 
high-performance  computing  center  collocated  at  Wright  Patterson  AFB.  For  typical  oxide- 
apertured  or  etched-post  VCSEL  structures,  the  RAM  required  was  typically  on  the  order 
of  1-3  GBytes.  Present  day  high  performance  desktop  stations,  however,  may  be  able  to 
implement  this  code  without  undue  strain.  Even  though  this  is  a  computationally  intensive 
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method,  this  model  has  the  distinct  advantage  that  it  calculates  the  total  optical  loss,  in¬ 
cluding  diffraction.  As  it  is  a  model  based  on  a  variatonal  solution  of  the  vector  Helmholtz 
equation,  its  results  not  only  allow  for  direct  calculation  of  lasing  mode  parameters,  but  also 
for  deeper  insight  into  the  underlying  physics  associated  with  VCSEL  design  parameters. 

In  the  sections  that  follow,  these  methods,  as  well  as  their  applications  (and  limitations!) 
to  VCSEL  design  will  be  detailed.  Furthermore,  a  section  that  describes  an  analytic  solution 
for  the  problem  of  graded  interfaces  between  two  slab  media  of  constant  refractive  indices 
will  also  be  presented.  This  solution  is  very  applicable  to  low-threshold  VCSEL  structures, 
as  much  of  the  wasted  power  in  driving  such  a  device  comes  from  trying  to  pass  current 
through  heterointerfaces.  The  use  of  material  grading  helps  to  overcome  much  of  this  waste, 
but  subsequent  optical  design  of  the  cavity  is  then  complicated  by  these  graded  layers. 

Weighted  Index  Method  with  Parasitics  (WIMP) 

One  goal  of  vertical  cavity  surface  emitting  laser  (VCSEL)  design  is  to  reduce  or  even 
eliminate  the  threshold  current.  In  relatively  large  devices,  the  threshold  current  is  pro¬ 
portional  to  the  volume  of  the  electrically-pumped  active  region.  Since  the  longitudinal 
dimension  of  this  volume  is  fixed  by  the  thickness  of  the  active  quantum  wells,  attention  has 
focused  on  transverse  current  confinement.  Both  etched-post  and  oxide- apertured  structures 
have  been  introduced  in  an  effort  to  restrict  the  current  path,  and  threshold  currents  have 
declined  accordingly  [4-6].  It  is  important  to  realize,  however,  that  these  lateral  patterning 
techniques  will  also  affect  the  transverse  optical  confinement.  As  the  transverse  dimensions 
shrink  to  the  order  of  the  lasing  wavelength — about  1  micron — two  optical  microcavity  ef¬ 
fects  emerge.  First,  the  spatial  profile  of  each  cavity  mode  changes  from  a  plane  wave  to 
a  true  three-dimensionally-confined  waveguide  mode.  This  change  alters  both  the  modal 
reflectivity  of  the  DBR  mirrors  and  the  modal  overlap  with  the  active  region  (the  transverse 
optical  confinement  factor).  Second,  the  energy  spacing  between  transverse  modes  increases 
sharply. 

Both  of  these  optical  microcavity  effects  could  significantly  influence  the  VCSEL  threshold 
current,  either  for  good  or  ill.  As  the  cavity  shrinks,  changes  in  the  modal  DBR  reflectivity 
and  the  transverse  optical  confinement  factor  will  modify  the  lasing  mode  threshold  gain. 
If  the  threshold  gain  increases  beyond  the  reach  of  the  quantum  wells,  the  structure  will 
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not  lase.  Conversely,  if  the  threshold  gain  remains  relatively  constant,  microcavity  VCSELs 
will  lase,  and  the  small  cavity  volume  will  result  in  a  low  threshold  current.  In  this  case, 
the  second  microcavity  effect  should  act  to  reduce  the  threshold  current  even  further:  as 
the  energy  spacing  between  modes  increases,  fewer  nonlasing  modes  exist  within  the  gain 
bandwidth,  and  therefore  fewer  carriers  recombine  via  nonlasing  modes.  These  submicron 
VCSELs  would  lase  at  some  of  the  lowest  currents  possible  in  semiconductor  lasers,  making 
them  attractive  for  system  applications.  They  would  also  exhibit  a  variety  of  quantum  elec- 
trodynamic  effects,  making  them  attractive  for  basic  laser  research  [7-9].  Since  fabricating 
submicron  VCSELs  is  difficult  and  expensive,  it  is  desirable  to  prove  the  effort  worthwhile 
by  first  calculating  the  lasing  thresholds  of  these  devices. 

Unfortunately,  the  bound  and  radiative  electromagnetic  modes  of  both  etched-post  and 
oxide- apertured  VCSEL  cavities  are  extremely  difficult  to  calculate.  Brute-force  numerical 
methods,  such  as  finite-difference  and  finite-element,  are  more  difficult  in  dielectric  structures 
than  in  metal-clad  waveguides  since  it  is  no  longer  possible  to  set  selected  field  components 
equal  to  zero  at  the  boundaries.  Absorbing  boundary  conditions  must  be  introduced  or  the 
structure  must  be  placed  inside  a  very  large,  perfectly-conducting  enclosure.  Furthermore,  in 
lasing  mode  calculations  we  must  solve  for  the  threshold  material  gain  in  the  active  region, 
meaning  that  the  dielectric  profile  of  the  structure  itself  is  one  of  the  unknowns.  These 
structural  unknowns  are  difficult  to  address  using  numerical  techniques. 

Often  it  is  better  to  address  problems  with  unknown  structural  parameters  through  an¬ 
alytic  calculations,  where  the  solutions  are  expanded  in  terms  of  an  infinite  set  of  special 
functions  and  the  expansion  coefficients  are  determined  by  boundary  conditions.  But  this 
technique  works  well  only  when  the  refractive  index  profile  separates  in  some  preferred  co¬ 
ordinate  system,  reducing  the  infinite  expansions  to  a  single  term.  If,  as  in  realistic  VCSEL 
structures,  the  refractive  index  profile  does  not  separate  in  any  coordinate  system,  then 
simple  single-term  special-function  solutions  to  the  governing  partial  differential  equations 
do  not  represent  exact  solutions  for  the  modes  and  analytic  methods  become  quite  cumber¬ 
some.  Despite  this  difficulty,  most  previous  analytic  calculations  have  introduced,  at  some 
point  in  the  treatment,  a  single  product  term  to  describe  a  particular  electromagnetic  field 
component  [10,  11].  It  is  vital  to  realize  that  this  is  equivalent  to  assuming  the  underlying 
differential  equation  separates. 

Since  separable  descriptions  facilitate  closed-form  expressions,  rapid  calculation,  and  com- 


parison  with  well- understood  “textbook”  problems,  there  is  considerable  motivation  to  im¬ 
prove  and  justify  them.  In  this  paper  we  generalize  the  weighted  index  method  (WIM) — a 
separable  approximation — to  compute  cavity  modes  in  cylindrically-symmetric  dielectric 
VCSEL  structures.  We  calculate  the  electric  and  magnetic  vector  potentials  and  use  these 
to  compute  the  resulting  fields.  We  also  show,  using  the  calculus  of  variations,  that  this 
technique  provides  the  best  separable  solution  to  the  scalar  Helmholtz  equation.  The  method 
allows  us  to  approximate  the  spatial  profile,  optical  confinement  factor,  resonant  frequency, 
and  threshold  material  gain  of  cavity  modes  in  both  oxide- apertured  and  etched-post  VC- 
SELs.  The  method  explicitly  considers  complex  media,  allowing  us  to  include  free  carrier 
losses.  We  start  by  summarizing  the  essential  vector  field  equations  needed  to  address  VC¬ 
SEL  modes.  Next  we  derive  the  WIM  and  outline  its  application  to  cylindrically-symmetric 
VCSEL  structures.  We  then  derive  weighted  boundary  conditions  needed  to  apply  the 
method  in  piecewise-constant  index  profiles.  Next,  we  outline  the  iterative  procedure  for 
solving  the  resulting  equations.  Finally  we  apply  the  method  to  oxide-apertured  and  etched- 
post  VCSELs,  computing  field  profiles,  optical  confinement  factors,  resonant  wavelengths, 
and  threshold  material  gains  for  several  cavity  modes.  A  summary  of  these  results  is  then 
presented. 

We  want  to  find  the  electric  and  magnetic  field  profiles,  the  resonant  wavelength 
(A),  and  the  threshold  material  gain  (gth)  for  each  cavity  mode  in  azimuthally-symmetric 
VCSEL  structures.  For  this  we  must  solve  a  vector-wave  equation  subject  to  appropriate 
boundary  conditions  at  each  interface.  Because  there  are  several  equivalent  electromagnetic 
descriptions  of  any  system,  we  can  write  wave  equations  for  the  electric  and  magnetic  fields, 
scalar  potentials,  or  vector  potentials.  The  most  powerful  and  convenient  method  for  this 
problem  is  to  solve  for  the  magnetic  and  electric  vector  potentials  and  use  them 
to  compute  the  fields.  The  steady-state,  time-harmonic  vector  potentials  A  and  F  satisfy 
the  three-dimensional  vector  Helmholtz  equation  (in  Gaussian  units) 


lu 


V  +-?e(p,z) 


A(p,  cj),  z) 
F(p,  <t>,  z) 


=  0. 


(1) 


Here  A  and  F  depend  on  time  as  elujt,  u o  =  27rc/A,  and  we  have  assumed  a  cylindrically- 
symmetric,  complex  dielectric  function  e.  Note  that 


\Je  =  N  =  n  +  in, 


(2) 
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where  N  is  the  (complex)  refractive  index;  material  gain  is  incorporated  by  taking  n  positive 
in  the  active  region.  We  assume  p  —  1  in  all  regions. 

The  power  of  the  vector  potential  approach  comes  from  the  fact  that  we  need  only  Fz  to 
generate  transverse  electric  (TE)  modes  and  only  Az  to  generate  transverse  magnetic  (TM) 
modes[12].  Since  an  arbitrary  electromagnetic  field  can  be  represented  as  a  superposition  of 
TE  and  TM  modes,  we  need  only  solve  for  the  two  unknown  scalar  functions  Az  and  Fz,  a 
dramatic  simplification  over  solving  Eq.  1  for  all  six  vector  components.  The  electric  and 
magnetic  fields  themselves  may  be  computed  directly  from  the  z  components  of  the  vector 
potentials  via  [13] 


E  =  —  - — V  x  V  x  (. zAz )  —  V  x  (zFz), 
H  =  V  x  (zAz) - V  x  V  x  (zFz), 


or,  more  explicitly, 


Fpip-j  z ) 

z) 

Ez{p ,  z) 
Hp(p,  z) 
H<p(p,z) 
Hz(p,z ) 


ic  d2  Id 

.  —  -£-Az(p,z)-~Fz(p,z), 

US  OpOZ  P  d(j) 


ic  f  d2  Id  m2 
cue  \  dp 2  ^  p  dp  p2 


Az(p,z), 


Id  icd2 

Ap' 


d  .  .  ic  d2 

op  (jjpp  ocpoz 


Ez  ( p ,  z) , 


IC 


d2  1  d 


m 


-  I  ~F~12  -~F) - 2  I  ^z(/p  z). 

cup  i  dp 2  pdp  pz 


(3) 

(4) 

(5) 

(6) 

(7) 

(8) 
(9) 

(10) 


Since  all  solutions  to  Eq  1  are  separable  in  the  azimuthal  coordinate  (f>,  depending  on  it 
as  the  differential  equation  for  Az  and  Fz  is  quite  simple.  By  expressing  A  and  F  in 
cylindrical  coordinates  and  inserting  the  appropriate  <f>  dependence  into  Eq  1,  we  have 


f  d2  Id  d2 

\  dp2  p  dp  dz2  ~*~ 


A  Z(p,z) 
Fz(p,  z) 


=  0. 


(11) 


In  separable  geometries,  the  two-dimensional,  azimuthally-symmetric  scalar  Helmholtz 
equation  (11)  may  be  solved  exactly  by  separation  of  variables,  yielding  the  potential  pro¬ 
files,  resonant  wavelength,  and  threshold  material  gain  for  each  cavity  mode.  Realistic 
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VCSEL  structures,  however,  are  not  separable,  greatly  complicating  the  solution  of  (11). 
An  exact  semianalytic  solution  could  be  obtained  by  expanding  Az  and  Fz  in  terms  of  the 
general  solutions  in  each  region,  then  matching  boundary  conditions  to  determine  the  (in¬ 
finite)  set  of  expansion  coefficients.  In  practice  this  technique  requires  considerable  care  to 
implement,  though  it  does  have  the  advantage  of  incorporating  non-separable  behavior  in 
the  solutions  [14] .  Below,  we  present  an  alternative  technique  to  generate  the  best  separable 
approximations  to  (11). 

Equation  (11)  represents  two  uncoupled  partial  differential  equations — one  each  for  Az 
and  Fz — which  are  quite  difficult  to  solve.  For  separable  geometries,  we  can  exactly  replace 
each  equation  in  (11)  with  two  independent  ordinary  differential  equations,  and  these  can 
be  solved  exactly.  For  non-separable  geometries,  we  approximate  the  solutions  to  (11).  In 
general,  there  are  two  possible  approximation  techniques.  The  most  common  approach  is 
to  maintain  the  exact  equations  (11)  and  construct  an  approximate  function  that  “almost” 
solves  them.  An  alternative  approach  is  to  replace  the  exact  equations  (11)  with  approximate 
equations,  and  solve  these  approximate  equation  exactly.  We  take  the  latter  approach,  and 
approximate  each  equation  in  (11)  with  two  coupled  ordinary  differential  equations.  We 
accomplish  this  by  extending  the  WIM — which  was  first  developed  to  calculate  waveguide 
modes  in  horizontal-cavity  ridge-waveguide  lasers  [15,  16] —  to  address  the  eigenmodes  of 
cylindrical  cavities.  This  technique  has  the  advantage  of  giving  the  best  separable  solution  to 
(11)  in  the  variational  sense,  and  allows  us  to  estimate  the  field  profile,  optical  confinement 
factor,  resonant  wavelength,  and  threshold  material  gain  of  each  cavity  mode.  Below  we 
derive  the  WIM  equations. 

Proceeding  as  if  separable  solutions  to  (11)  exist,  we  take 


Az(p,  z)  =  P(p)Q(z), 
Fz(p,z)  =  R(p)S(z). 

Substituting  either  of  these  into  (11)  gives 


(12) 


C "(pm  +  -  A)  =  0,  (13) 

p  V c  p  J 

where  (  =  P  or  R  and  £  =  Q  or  S.  For  each  potential  in  (12),  we  can  separate  the 
resulting  equation  (13)  by  integrating  it  against  (*(p)  or  This  procedure  yields  the 
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axial  equation 


rw  +  few)bw  =  ° 


(14) 


and  the  radial  equation 


C"M  +  k'M  + 


21 


few)2  -  y 


C  (p)  =  o 


(15) 


for  each  potential.  These  axial  and  radial  equations  are  coupled  by  the  weighted  axial  and 
radial  propagation  constants,  given  respectively  by 


few)  = 


2_  W2/c2)(C|£W^)|C)4«IC")  +  (C|p-1IC')-m2{Cli>-2IC) 


and 


(CIO 


(ha  ,  w2_02/c2)(C|s(pfe)|0  +  (Cin 


Here  (  |  )  denotes  an  inner-product  over  z  or  p,  respectively  defined  by 

/OO 

A*(z)B(z)  dz 


(16) 


(17) 


(18) 


and 

POO 

(A(p)\B(p)}=  A*(p)B(p)pdp.  (19) 

do 

Since  the  weighted  propagation  constants  depend  on  whether  we  solve  for  Az  or  Fz,  we 
have  introduced  an  additional  superscript  a  =  TE  or  TM  to  distinguish  between  TE  modes 
resulting  from  Fz  (involving  averages  over  (  =  R  and  £  =  S)  and  TM  modes  resulting  from 
Az  (involving  averages  over  (  =  P  and  £  =  Q).  Coupling  occurs  only  between  the  radial  and 
axial  equations  for  a  given  vector  potential  Az  or  Fz:  the  two  potentials  remain  uncoupled 
in  (14)  and  (15).  But  Az  and  Fz  will  be  coupled  later  by  boundary  conditions  when  we  solve 
for  hybrid  modes. 

We  must  now  find  the  solutions  of  (14)  and  (15)  in  piecewise-constant  refractive  index 
profiles,  paying  particular  attention  to  the  interfacial  boundary  conditions.  A  sample  struc¬ 
ture  is  shown  in  Fig.  2.  Since  we  will  work  exclusively  with  piecewise-constant  geometries, 
we  simplify  our  notation  by  taking  e(p,  z)  —?  where  i  and  j  index  the  radial  and  axial 
regions,  respectively.  Thus  we  also  have  k^(p)  — >  and  /3“//  (z)  /?“,  and  (14)  and  (15) 
reduce,  respectively,  to  the  one- dimensional  Helmholtz  equation  and  Bessel’s  equation. 

To  solve  these  equations  subject  to  the  refractive  index  profile  c,tJ  of  etched-post  or 
oxide-apertured  VCSEL  structures,  we  must  supplement  (14)  and  (15)  with  an  appropriate 
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FIG.  2:  Illustration  of  piecewise  constant  permittivity  notation  for  an  axially  symmetric  VCSEL. 


set  of  interface  and  endpoint  boundary  conditions.  The  interfacial  boundary  conditions 
are  the  usual  continuity  requirements  on  the  normal  and  tangential  components  of  various 
electromagnetic  fields.  Therefore  we  cannot  directly  enforce  boundary  conditions  on  P,  Q, 
R,  and  S,  but  must  perform  the  intermediate  step  of  computing  the  electric  and  magnetic 
fields  via  (5)-(10).  Furthermore,  since  the  underlying  partial  differential  equations  (11)  do 
not  separate,  these  boundary  conditions  cannot  be  satisfied  at  all  points  on  the  boundary 
surfaces-if  they  could,  (11)  would  be  separable. 

In  preparation  for  generating  approximate  boundary  conditions,  we  rewrite  the  weighted 
index  formulas  (16)  and  (17).  By  integrating  (15)  against  ( *(p)  we  have 

(CIC")  +  (C|p-1IC/)-m2(C|p-2IC)  =  -(CI(^)2IC),  (20) 


allowing  us  to  express  (16)  as 


foaP  (^2/c2)(  C  I  £ij  |  C  )  -  (  C  I  (K)2  K  }  _  ^  ,  a  ,  / 1  a  \  2 
\Pj  )  -  -  \K  )  ’ 


where 


^'-(CkylC)  md  (ka)=  /<CI  (fcf)UC) 
{£u)-  <  C I  c  >  {k)-\  (CIO 


Similarly,  by  integrating  (14)  against  £*(z)  we  find 


<£ir>  =  -ui(/3?m, 


allowing  us  to  express  (17)  as 

ft")2  = 


2  (^y2){?k„u>-«;i(y)2io 


uio 


(21) 


(22) 


(23) 


(24) 
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where 


>  = 


_  ( { I  £»,j  I  £ ) 
(£10 


and  (/?“)  = 


(£1  (/5“)2IO 

(eio 


(25) 


The  compact  expressions  (21)  and  (24)  allow  us  to  compute  weighted  variables  without  using 
the  derivatives  of  (  and  £. 

We  specialize  these  expressions  to  VCSEL  lasing  modes  by  truncating  the  inner  product 
over  z,  replacing  (18)  with 


/‘2-max 

( A(z)  |  B(z) )  =  /  A*(z)B(z)  dz,  (26) 

**  2-min 

where  zm{n  and  zmax  denote  the  lower  and  upper  VCSEL  boundary  planes.  This  truncation 
is  necessary  to  force  (24)  to  converge  and  reflects  the  assumption  that  most  of  the  energy  is 
contained  inside  the  VCSEL  cavity.  In  contrast,  we  force  the  radial  wavefunctions  to  decay 
evanescently  to  zero — to  find  guided  modes — and  the  inner  product  defined  by  (19)  presents 
no  difficulty. 

We  now  present  boundary  conditions  and  solutions  for  the  axial  and  radial  equations 
(14)  and  (15).  Recall  that  the  unknowns  in  this  formalism  are,  for  each  mode,  the  functions 
P(p),  Q(z),  R{p), and  S(z),  the  resonant  frequency  cu,  and  the  threshold  material  gain  in  the 
active  region  gth  =  ^Kactive/ A. 

We  solve  the  axial  equation  (14)  in  piecewise  constant  geometries,  such  as  in  Fig.  3.  The 
general  solutions  of  (14)  are  given  in  each  axial  region  [z3 ,  z3+ 1]  by 


Q3(z) 


=  aJMe‘ ^ 


+  b 


TM„-ipJMz 


(27) 


Sj{z)  =  aJEei0f ’’  +  bJBe-igr*, 

where  we  explicitly  denote  both  the  TE  and  TM  solutions  for  clarity.  Using  the  iterative  so¬ 
lution  procedure  described  below,  we  compute  /3jE  and  /3jM  from  (21);  assume  for  now  that 
they  are  known  constants.  These  solutions  must  be  joined  at  each  interface  Zj  by  matching 
the  tangential  electric  and  magnetic  fields.  Inserting  (12)  into  (5)— (10),  we  compute  these 
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FIG.  3:  Illustration  of  weighted  permittivity  profile  for  the  WIM  axial  solution. 


tangential  field  components  as 

i  n  i  rri 

Ep(p,  z)  = - P\p)Q\z)  -  — R(p)S(z),  (28) 

ujs  p 

run 

E+(p,  z)  =  — P(p)Q'(z )  +  R'(p)S(z),  (29) 

usp 

i  m  i  c 

Hp(p,  z)  =  — P{p)Q{z)  -  —R\p)S\z),  (30) 

p  Id  pi 

cm 

z)  =  - P'(p)Q(z )  + - R{p)S\z).  (31) 

(jjpp 

We  consider,  in  turn,  two  distinct  cases:  m  =  0  and  m^O. 

When  m  =  0  we  can  match  boundary  conditions  with  Q  =  0  (pure  TE  modes)  or 
S  =  0  (pure  TM  modes).  For  these  modes,  it  is  sufficient  to  force  just  two  tangential 
field  components  to  be  continuous:  demanding  continuity  of  the  other  two  components  gives 
redundant  conditions.  For  TE  modes  we  require  E,p  and  Hp  to  be  continuous;  for  TM  modes 
we  require  Ep  and  H0  to  be  continuous.  As  we  will  see  below,  for  pure  TE  and  TM  modes 
the  fields  themselves,  and  not  just  the  vector  potentials,  are  separable.  This  makes  it  easy 
to  generate  weighted  boundary  conditions  for  these  modes. 

For  pure  TE  modes,  we  have 


E*{p,z) 

=  K(p)S(z), 

(32) 

Hp(p,  z) 

=  ~—R'(p)S\z). 

(jJ  P 

(33) 
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Since  we  assume  p  =  1  in  all  regions,  both  E,p  and  Hp  depend  on  z  only  via  S(z)  and  S'(z), 
respectively.  Therefore  we  can  force  both  tangential  fields  to  be  continuous  by  setting  S(z) 
and  S'(z)  continuous  across  each  interface. 

For  pure  TM  modes,  on  the  other  hand,  we  have 


1  c 

Ep(p,z )  = - P'(p)Q\z), 

UJS 

(34) 

H^Pi  z)  =  P'(p)Q(z). 

(35) 

Again,  we  can  make  continuous  by  forcing  Q(z)  to  be  continuous  across  each  interface. 
However,  Ep  depends  on  2  through  both  Q’(z)  and  e.  Since  e  is  not  a  separable  function 
of  2  and  p,  we  must  weight  the  relative  permittivity  to  obtain  average  boundary  conditions 
holding  for  all  p.  Therefore  we  require  Q\z)/(epM)  to  be  continuous  at  each  interface, 
where  (e™)  has  already  been  defined  in  each  region  by  (25).  (We  could  also  have  generated 
average  boundary  conditions  using  (^),  but  this  approach  gave  inferior  results.) 


Inserting  the  functional  forms  (27)  and  applying  the  continuity  conditions  for  either  a  = 
TE  or  a  =  TM  modes,  we  relate  a",  6"  to  a“+1,  6“+1  at  each  axial  boundary  z  =  Zj  through 
the  transfer  matrices 


The  composite  transfer  matrix  [Ta]  for  the  whole  system  is  formed  by  cascading  the  indi¬ 
vidual  transfer  matrices,  giving 


[T“]  =  [H-'KPr'KHi? 


where  N  is  the  number  of  axial  regions  (including  substrate  and  air)  in  the  problem  geometry. 
Thus  we  relate  the  unknown  coefficients  in  the  j  =  1  region  (substrate)  to  the  coefficients 
in  the  j  —  N  region  (air)  via 
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Since  we  are  searching  for  axially-emitting  (lasing)  modes,  we  permit  only  outgoing  radiation 
by  setting 

b«  =  aaN  =  0.  (41) 

Finally,  substituting  (41)  into  (40)  we  obtain  the  axial  threshold  condition 

^22 iP  i  ^ active )  0.  (42) 


Setting  the  real  and  imaginary  parts  of  t%2  equal  to  zero  gives  two  independent  equations 
that  we  solve  to  obtain  the  modal  frequency  u  =  2irc/\  and  threshold  material  gain  gth  = 
YnKactlve/ A  for  pure  TE  and  TM  modes.  The  expansion  coefficients  a"  and  6"  for  each  region 
are  found  by  back  substitution  through  (40)  and  (36). 

In  order  to  generate  sensible  boundary  conditions  when  m  ^  0,  we  must  construct  hybrid 
modes  in  which  both  S  and  Q  nonzero.  In  this  case  none  of  the  tangential  fields  are  separable, 
since  each  has  both  TE  and  TM  parts.  Each  TE  and  TM  part  is,  in  turn,  a  sum  of  cylindrical 
wave  terms  like  eirn^Jm(kp)ell3z  (P  and  R  will  turn  out  to  be  Bessel  functions).  An  exact 
solution  would  require  us  to  include  a  superposition  of  cylindrical  waves  involving  all  values 
of  k,  and  the  boundary  conditions  would  couple  all  terms  at  each  interface.  Fortunately,  the 
coupling  between  cylindrical  waves  with  different  k  is  small  enough  to  ignore  [14,  17,  18]. 
Therefore  we  simultaneously  and  independently  enforce  the  continuity  of  the  TE  and  TM 
mode  components  of  each  tangential  field,  and  our  hybrid  mode  boundary  conditions  become 
the  same  as  those  for  pure  TE  and  TM  modes. 

We  solve  the  radial  equation  (15)  in  piecewise  constant  geometries,  such  as  in  Fig.  4. 
The  general  solutions  of  (15)  are  given  in  each  radial  region  [pi}  pi+i]  by 


cfMJm(kfMp )  +  dJMYm(kfMp)  i  ±  M 
clMKm(ikTMMp)  +  dTMMIm(ikTMMp )  i  =  M  ’ 

CJE  Jm(kfE  p)  +  dfEYm(k'[E  p)  i^M 
c^fKm(ikJfp)  +  dlflm(ik^Ep)  i  =  M 


(43) 


where  Jm  and  Ym  are  m-th  order  Bessel  functions  of  the  first  and  second  kind,  Im  and  Km 
are  modified  m-th  order  Bessel  functions  of  the  first  and  second  kind,  and  i  =  1,2 , ,M 
indexes  the  inner  to  outer  radial  regions.  Using  the  iterative  solution  procedure  described 
later,  we  compute  kfE  and  kfM  from  (24);  assume  for  now  that  they  are  known  constants. 
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FIG.  4:  Illustration  of  weighted  permittivity  profile  for  the  WIM  radial  solution. 

These  solutions  must  be  joined  at  each  interface  pj  by  matching  tangential  electric  and 
magnetic  fields.  Inserting  (12)  into  (5)-(10)  we  compute  these  tangential  held  components 
as 


E^{p,z)  = 


P(p)Q(z), 

uoe 

(44) 

run 

p(p)Q'{Z)  +  R'(p)S{z), 

ujsp 

(45) 

(kfE)2  R(p)S(z), 

(jJ  /i 

(46) 

cun 

-P\p)Q{z)  + - R(p)S\z). 

(jjpp 

(47) 

We  consider,  in  turn,  the  distinct  cases  m  =  0  and  m^O. 

As  in  the  axial  problem,  when  m  =  0  we  can  match  boundary  conditions  with  P  =  0 
(pure  TE  modes)  or  R  =  0  (pure  TM  modes),  and  it  is  sufficient  to  force  just  two  tangential 
held  components  to  be  continuous.  For  TE  modes  we  require  E,p  and  Hz  to  be  continuous; 
for  TM  modes  we  require  Ez  and  H4>  to  be  continuous.  Following  the  same  arguments  as  in 
the  axial  problem,  we  construct  weighted  boundary  conditions  to  require  the  continuity  of 


{kfE)2  R(p)  and  R'(p )  (for  TE  modes), 
(k™)2 

/  ~tW\  R(p )  and  P\p)  (f°r  TM  modes), 

\£q  ) 


(48) 

(49) 
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at  each  radial  interface  p  =  p*.  To  generate  radial  boundary  conditions  independent  of  z, 
we  have  replaced  £hj  by  the  weighted  permittivity  (c™)  defined  in  (  25). 

Inserting  the  functional  forms  (43)  and  applying  the  continuity  conditions  for  either  a  = 
TE  or  TM  modes,  we  relate  c“,  d"  to  c“+1,  df+1  at  each  radial  boundary  p  =  pi  through  the 
transfer  matrices 


ra 
^ i 

=  b? 

ci+ 1 

l - 

a- 

rfa 

ai+ 1 

where 


A 


TE  — 


TE  — 


B 


aTm  = 


BjE  = 


A? 


(kfE)2  Jm{kJEpi )  (kJE)2  Ym{kJE pi) 
kJEJ'm(kJEPi)  kJEY}n{kJE  pi) 

{kJ+iY  Emi'lfkTftpi)  (klfi)2  UmilkT+iPi) 
AkJ+i  (7  kEE{  pi )  lkEEx  H'm  (7  kEE{  pt ) 


(50) 


fc™)8  7  {*T“Y 


■ ) 


Ym(kIMPi) 


kfM  J'm{kJM  Pi)  k™Y^k™Pi) 


h™  , 

SfrLn UlklXp.) 

\S+i,g/  \£z+i,q/ 

7  k™  Bm  (7^™  Pi )  7  k™  Wm  (7^™  Pi ) 


k™' 
l  Ki+ 1  / 


(51) 

(52) 

(53) 

(54) 


Here  we  have  defined 


7  = 


1  inner  regions 
outer  region 


Jm  inner  regions 
Km  outer  region 


n. 


As  in  (39),  we  form  a  composite  system  matrix 


\ua]  =  Kl'Msa  •  •  •  [Ah- ,rK-,] 


(56) 


relating  the  unknown  coefficients  in  the  i  —  1  region  (core)  to  the  coefficients  in  the  i  =  M 
region  (cladding)  via 

(57) 

Since  we  are  searching  for  longitudinally-propagating,  laterally-confined  VCSEL  modes  we 
force  regularity  at  the  origin  and  exponential  decay  as  p  — »  00  by  setting 


nOi 

C1 

un 

<2 

1 

O 

l 

a- 

1 

£ 

top 

U22 

1 

d“  =  dB  =  0. 


(58) 
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Finally,  substituting  (58)  into  (57)  we  obtain  the  radial  threshold  condition 

^21  ^ active )  0.  (59) 

In  order  to  generate  sensible  boundary  conditions  when  m  ^  0,  we  must  construct  hybrid 
modes  in  which  both  R  and  P  are  nonzero.  Again,  none  of  the  tangential  field  components 
are  separable  and  it  is  impossible  to  match  them  for  all  z  at  a  radial  interface.  Further¬ 
more,  if  we  mirror  our  axial  treatment  and  independently  force  the  TE  and  TM  components 
to  be  continuous  we  generate  an  inconsistent  boundary  condition.  For  example,  indepen¬ 
dently  forcing  Ez  and  the  TM  part  of  E,p  to  be  continuous  requires  ( kfM )2  P(p)/ (e™)  and 
P(p)/ (sqM)  be  continuous  at  each  interface.  But  this  forces  kfM  itself  to  be  continuous, 
which  is  absurd.  We  cannot  hope  to  work  with  TE  and  TM  modes  independently.  Therefore 
we  approximate  our  problem  with  an  equivalent  cylindrical  dielectric  waveguide  problem, 
which  admits  analytic  solutions  for  hybrid  modes.  These  approximations  are  used  only  to 
generate  boundary  conditions,  not  as  a  substitution  for  the  actual  radial  and  axial  solutions 
in  each  region. 

Cylindrical  waveguide  modes  depend  on  z  as  el^z  or  e~i/3z.  Our  fields,  in  contrast,  depend 
on  z  via  (27).  There  are  two  differences  we  must  overcome.  First,  cylindrical  waveguide 
modes  are  characterized  by  a  single  axial  propagation  constant  /3,  whereas  our  fields  have  a 
different  axial  propagation  constant  /?“  for  each  axial  region  and  polarization.  We  can  easily 
remedy  this  by  replacing  /?“  with  (/?“),  as  defined  in  (25).  Second,  cylindrical  waveguide 
modes  depend  on  z  as  either  etzlz  or  e~l,iz .  But  even  after  replacing  /?“  with  (fP)  our  fields 
have  a  different  linear  combination  of  el^a')Z  and  e~l^a',z  in  each  region,  depending  on  the 
relative  values  of  a"  and  6".  Therefore,  for  the  purpose  of  constructing  boundary  conditions, 
we  assume  that  the  fields  approximate  a  “pure”  standing  wave  in  the  axial  direction,  with 

aJE  =  bEE  and  aJM  =  —bJM.  (60) 

We  assume  further  that  aJE  ~  aJM ,  and  approximate  the  ^-dependence  of  our  fields  as 

Q(z)  =  sin  ((P™)z), 

(61) 

S(z)  =  cos  ((/3TE)z). 
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Substituting  (61)  into  (45)  and  (47),  we  approximate  the  0  field  components  as 


Mp,z)  =  -  TM  P(p)Q'(z)  +  R'(p)S(z), 

)P 

«  I  P(P)  +  R'(P)  I  cos((/3>,)  (62) 

and 

~  ~  cm  ~ 

#*(**)  =  -P\p)Q{z)  +  — i2(p)5'(z), 

(jJ  P 

^  — -#(/?))  sin  ((/?)z) .  (63) 

l  cup  J 

Here  we  have  replaced  c,tJ  by  the  appropriate  weighted  value,  and  have  selectively  assumed 
(/?TE)  «  {0™)  ~  (/?)  to  factor  out  the  z  dependence.  Through  a-posteriori  comparison 
with  our  calculated  results,  we  find  all  these  assumptions  well  justified. 

These  approximate  expressions  for  E(p  and  are  separable,  as  are  expressions  (44)  and 
(46)  for  Ez  and  Hz.  Therefore  we  can  immediately  construct  suitable  boundary  conditions 
for  radial  hybrid  modes  by  requiring  the  continuity  of 

( u TM)2 

^myp^  (64) 

(kJE)2R(p),  (65) 

CTYll  3^^') 

-^M^P(P)+R'(P)’  (66) 

cm  ( 3^E  \ 

and  P\p)  +  1  R(p).  (67) 

LJPOP 


By  inserting  the  functional  forms  (43)  and  forcing  all  four  of  the  above  combinations  to  be 
continuous,  we  link  the  unknown  coefficients  at  each  radial  interface  p  =  p*  through  the  4x4 


transfer  matrix 


’  „TM  " 

"  CTM  " 

Ci 

Ci+1 

dJM 

=  B% 

JTM 
ai+ 1 

nTE 

„TE 

Ci 

Ci+ 1 

AE . 

dTE 

L  ai+ 1  J 

(68) 
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where 


A* 


O _ L  J  (k™n  ) 

I jt M \  Hi) 

\£i,Q  / 

0 


(*? 


0 


cm(P™)  j  /  uTM n  \  cm(0™)v  (uTMn\  h,TE  V  (uTEn\  kTEV  (hTEn\ 
uj(£Tq)p  'll  '  u(e™)p  Pi)  hi  Jm\Ei  Pi)  A  1  m\Ei  Pi) 


Bi  = 


(kfE)  Jm(kfEPi)  ( kIE )  Ym(kfEPi) 
kJE  UkJE Pi)  kJEYJn(kJE  Pi) 

UTM  77  (uTMn\  h.TMV'  (k™  n\  cm(/AE)  7  (kTEn\  7  (kTEn\ 

Hi)  1  m\^i  Hi)  up  Jm\^i  Hi)  up  Jm\^i  Hi) 

(fc»™)2-  (  kTM  \  (fc,™)2-Q  (  kTM  \ 


(69) 


2  +  1  ' 


i+1,  Q 


0 


0  0 

i  2  i — i  /  i  T'TP.  \  /  j  TR  \  2 


(*w)  Sm(7^p*)  (^)  nm(7^Pi) 
+^^-sm(7fc™P7)  +^+nm(7fc™pi)  7fc™s^(7/c^pi)  7fc™n^(7fc™pi) 


cm(/3™) 


_  ^k™E,m('yk™  pp  ikT^U'^kHf  pP 


cm((3TE)  ^ 


u;pi 


cm(/3TE) 


LOPi 


nm(7*mPi) 

(70) 


These  matrices  are  mathematically  equivalent  to  those  for  a  cylindrical  dielectric  waveguide. 
Cascading  the  interface  transfer  matrices,  we  again  derive  a  composite  system  transfer  matrix 


At  TM  " 
C1 

un 

U\ 2 

«13 

Uli 

1 

C7> 

d™ 

«21 

«22 

«23 

U24 

JEM 

aM 

rTE 

C1 

W.31 

W32 

m33 

U34 

rTE 

CM 

_d\E  _ 

_  ^11 

«42 

W43 

ii44 

dTE 

L  aM  J 

(71) 


relating  the  innermost  and  outermost  radial  coefficients.  The  endpoint  boundary  conditions 
(58)  remain  valid.  Applying  them  to  (71)  and  demanding  nontrivial  solutions  gives  the 
hybrid  threshold  condition 

W21  {id ,  At  active  )  ^23  (^ ;  & active  ) 

(i  ll  (uJ .  H  active)  2/43  Inactive ) 

which  we  solve  in  the  complex  plane  for  c o  and  nactlve . 

We  compute  the  longitudinal  and  transverse  mode  spectrum  by  self-consistently  solving 
the  radial  and  axial  problems.  The  modes  are  specified  by  the  longitudinal  mode  number, 
the  transverse  mode  number,  and  the  azimuthal  mode  number  to.  Different  longitudinal 
modes  correspond  to  successive  roots  of  the  axial  threshold  condition  (42),  while  different 
transverse  modes  correspond  to  successive  roots  of  the  radial  threshold  conditions  (59)  (for 
TE/TM  modes)  or  (72)  (for  hybrid  modes).  The  energy  spacing  between  longitudinal  modes 
is  much  greater  than  that  between  transverse  modes.  If  we  let  n  denote  a  generalized  mode 


=  0, 


(72) 
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index  corresponding  to  a  particular  TE,  TM,  or  hybrid  mode,  then  the  explicit  unknowns  for 
each  mode  are  the  vector  potential  functions  Pn(p),  Qn(z)  and/or  Rn(p),  Sn(z ),  the  optical 
mode  frequency  u>n,  and  threshold  values  for  (, inactive) n ■  We  iteratively  solve  the  axial  and 
radial  problems  as  follows. 

First  we  focus  attention  on  a  particular  family  of  modes  by  fixing  m:  for  m  =  0  we  can 
compute  TE  or  TM  modes,  while  for  m  >  0  only  hybrid  modes  are  possible.  Then  we  solve 


the  axial  problem,  taking  the  effective  indices  equal  to  the  corresponding  values  for  er  in 
the  innermost  radial  region:  this  solution  corresponds  to  the  standard  plane- wave  calculation 
appropriate  for  large-area  devices.  The  ordered  roots  of  the  axial  threshold  condition  (42) 


yield  initial  approximations  Qno{z)  and/or  Sno(z),  ujn0xml\  and  (tractive  )  f°r  the  mode  n0 
of  interest.  For  VCSEL  lasing  mode  calculations  we  are  only  interested  in  the  first  root 
of  (42),  corresponding  to  the  fundamental  longitudinal  mode.  We  next  compute  the  kf 
by  inserting  Qno(z )  and/or  Sno(z)  into  (24),  giving  us  enough  information  to  address  the 
radial  problem.  Depending  on  whether  we  are  solving  for  TE,  TM  or  hybrid  modes,  we  then 
compute  the  roots  of  the  radial  threshold  condition  (59)  or  (72),  yielding  approximations  for 


axial ) 


Pno(p )  and/or  Rno(p),  u>n“dml\  and  (  )  :  successive  roots  correspond  to  progressively 

higher-order  transverse  modes.  Then  we  alternate  between  solving  the  axial  and  radial 
problems,  always  updating  /?“  and  kf  by  inserting  the  most  recent  wavefunctions  into  (21) 
or  (24)  as  appropriate.  In  this  way  we  generate  a  self-consistent  solution  to  the  coupled 
WIM  equations  (14)  and  (15),  terminating  when 


_( radial ) 


(axial)  ( radial ) 

o  5 

f(axial)\  _  (  (radial)\ 

\  ™ active  I  \  ™ active  I 

\  'no  \  /no 

to  within  prescribed  tolerance.  The  procedure  converges  quite  rapidly,  allowing  us  to  solve 
for  a  large  number  of  cavity  modes. 

Finally,  we  note  that  the  character  of  the  modes  found  depends  on  both  the  differential 
equations  and  on  the  endpoint  boundary  conditions  enforced.  The  original  application  of  the 
WIM  to  rectangular  waveguide  geometries  assumed  propagating  behavior  in  the  z  direction 
and  evanescent  decay  in  the  x  and  y  directions.  These  endpoint  conditions  resulted  in  a  stan¬ 
dard  eigenvalue  problem,  with  the  longitudinal  propagation  constant  (3  being  the  eigenvalue. 
In  this  case,  the  Rayleigh- Ritz  variational  principle  asserts  that  the  resulting  approximation 
of  f3  will  be  more  accurate  than  the  wavefunctions  themselves.  In  our  application,  we  have 
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enforced  evanescent  decay  in  the  radial  direction  and  have  permitted  only  outward  propa¬ 
gating  waves  in  the  z  direction.  This  constrains  the  mode  frequency  u  and  the  threshold 
material  gain  n  ,  as  opposed  to  the  propagation  constant  / 3 ,  and  these  unknown  parameters 
no  longer  appear  as  eigenvalues.  Therefore  the  Rayleigh- Ritz  principle  has  nothing  to  say 
about  the  relative  accuracy  of  u  and  k.  Nevertheless,  general  variational  principles  dictate 
that  we  have  found  the  best  separable  solution,  and  therefore  we  expect  the  resulting  values 
for  c o  and  k  to  be  reasonably  accurate.  Note  that  we  could  also  have  solved  for  radial  lasing 
modes  by  requiring  evanescent  decay  in  the  axial  direction  and  outward  radial  propagation, 
though  these  would  be  difficult  to  approximate  because  of  the  infinite  radial  inner  product. 
In  fact,  as  the  cavity  radius  shrinks  we  do  see  evidence  of  radial  propagation,  as  discussed 
in  the  next  section. 

Using  the  methods  described  above,  we  calculate  lasing  modes  in  etched-post  [19,  20] 
and  oxide-apertured  [5,  21]  devices  fabricated  from  a  A  =980  nm,  1.5A-cavity  VCSEL.  The 
VCSEL  reflectors  consist  of  a  17.5  period  p-type  GaAs/Alo^Gao.osAs  top  Distributed  Bragg 
Reflector  (DBR)  for  the  oxide-apertured  structure,  and  a  4  period  GaAs/ARO  top  DBR  with 
a  A/4  p+  GaAs  contact  layer  for  the  post  geometry.  The  bottom  reflector  in  both  struc¬ 
tures  is  a  22  period  n-type  GaAs/Alo^Gao.osAs  DBR.  The  1.5A-cavity  is  step-tapered  with 
(intrinsic)  layers  of  Alo.9sGao.02As,  Alo.65Gao.35As,  Al0.30Ga0.70As,  and  GaAs,  culminating  in 
a  single  In0.2Ga0.sAs  quantum  well.  Both  structures  are  grown  on  a  GaAs  substrate.  The 
etched-post  structure,  illustrated  in  Fig.  5,  is  formed  by  etching  the  top  GaAs/AlAs  DBR 
down  to  the  A/4  GaAs  contact  layer,  then  oxidizing  the  AlAs  layers.  The  oxide-apertured 
structure,  illustrated  in  Fig.  6,  is  formed  by  oxidizing  the  A/4  Al.9sGa.02 As  layers  in  the 
cavity.  We  model  each  region  as  a  cylindrically-symmetric  layer  of  constant  refractive  index, 
assuming  the  material  parameters  in  Table  I.  Free  carrier  losses  are  incorporated  by  taking 
Kitj  negative.  Material  gain  is  incorporated  by  taking  Kactlve  positive  in  the  active  quantum 
well  region. 

The  lowest  frequency,  or  fundamental,  VCSEL  lasing  mode  is  analogous  to  the  HEMn 
(or  HEn)  hybrid  waveguide  mode.  Here  the  first  and  second  subscripts  refer  respectively 
to  the  azimuthal  (to)  and  radial  mode  numbers.  The  HEMn  mode  is  the  most  plane-wave 
like  of  all  the  propagating  bound  modes,  despite  containing  both  Ez  and  Hz  [22,  23].  It 
is  also  the  only  waveguide  mode  having  a  radial  intensity  distribution  with  a  maximum  at 
the  center.  Following  this  terminology,  we  refer  to  the  next  higher-order  VCSEL  modes  as 
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FIG.  5:  Quasi-3D  plot  of  the  etched-post  YCSEL  index  profile. 

HEM2i,  TE0i,  and  TM0i  modes:  these  modes  make  up  the  degenerate  LPn  mode  under  the 
“linear  polarization”  approximation  [23] .  All  three  of  these  modes  feature  a  radial  intensity 
profile  with  a  null  at  the  center. 

In  Figs.  7  and  8  we  plot  the  longitudinal  held  profile  for  sample  etched-post  and  oxide- 
apertured  VCSELs.  For  the  oxide-apertured  structure,  the  top  DBR  is  very  long  and  the 
fields  penetrate  deeply  into  the  top  DBR.  Therefore  the  large  index  difference  between 
the  semiconductor  cavity  and  surrounding  oxide  is  heavily  weighted  in  (25),  giving  a  large 
discontinuity  in  (e^)  between  the  inner  and  outer  radial  regions.  As  a  result,  the  fields 
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FIG.  6:  Quasi-3D  plot  of  the  oxide-apertured  VCSEL  index  profile. 


within  the  oxide-aperture  VCSEL  are  tightly  confined  to  the  core.  In  contrast,  for  the 
etched-post  structure  the  large  index  contrast  between  top  DBR  layers  allows  very  little 
field  penetration,  resulting  in  a  smaller  effective  index  difference  and  a  correspondingly  less 
confined  field.  In  both  structures,  the  fundamental  mode  is  more  confined  than  higher-order 
modes,  and  smaller  VCSELs  exhibit  less  confinement  than  larger  ones.  These  effects  can 
also  be  seen  in  the  three-dimensional  field  intensity,  given  by  [24] 


2 

2" 

£ 

E 

+  /X 

H 

(73) 


which  we  plot  in  Figs.  9-10.  The  radial  discontinuity  in  these  figures  is  a  measure  of  the 
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TABLE  I:  Material  Parameters  used  for  WIM  Simulations 


Material 

Index  Doping  (1018  cm  3) 

Loss  (cm  1) 

ino.2Gao.8As 

3.5691 

none 

N/A 

GaAs 

3.5256 

none 

N/A 

Alo.3Gao.7As 

3.3622 

none 

N/A 

Alo.65Gao.35As 

3.1637 

none 

N/A 

Al0.98Ga0.02As 

2.9713 

none 

N/A 

AlxO 

1.55 

none 

N/A 

n-GaAs 

3.5256 

2 

10 

n-Alo.92Gao.08As  3.0067 

2 

10 

p-GaAs 

3.5256 

1 

11.5 

p-Alo.92Gao.08As  3.0067 

1 

11.5 

p+-GaAs 

3.5256 

5 

57.5 

error  in  our  solution.  This  error  stems  from  our  separable  approximation. 

The  transverse  confinement  factor  T4  is  usually  defined  as 

f  ..  \E\2  ds 

p  J  active  1  1 _ 

J\E\2ds  ’ 

where  the  integral  in  the  numerator  is  over  the  transverse  extent  of  the  active  region  and  the 
integral  in  the  denominator  is  over  the  entire  transverse  extent  of  the  field.  Figure  11  shows 
estimates  of  the  transverse  confinement  factor  for  the  first  two  modes  of  the  etched-post  and 
oxide-apertured  VCSELs  as  a  function  of  cavity  radius.  The  estimates  are  generated  from 
(74),  using  E  =  E(p)  — ■>  E^(p)  for  TE  modes  (see  (62))  and  E  — >  H(p{p)  for  TM  modes  (see 
(63));  for  hybrid  modes  we  use  both  E^{p)  and  H,p(p)  and  average  the  results.  We  use  the 
effective  <p  components  since  they  are  already  averaged  over  z  and  are  representative  of  the 
total  field  intensity  profile.  As  the  cavity  radius  decreases,  more  of  the  field  intensity  leaks  out 
of  the  active  region  and  the  confinement  factors  drop  monotonically.  This  behavior  becomes 
more  pronounced  for  higher-order  modes,  as  illustrated  in  Fig.  12.  Our  oxide-apertured  de¬ 
vices  confine  the  optical  mode  more  strongly  to  the  active  region  and  have  higher  transverse 
confinement  factors  than  our  etched-post  structures.  But  the  larger  rate  of  change  in  confine- 


27 


Position  (jam) 


FIG.  7:  Index  and  standing  intensity  profile  along  the  axial  direction  for  a  1.4  pm  radius  etched- 
post  VCSEL. 

ment  factor  for  the  etched-post  VCSEL,  as  shown  in  11,  yields  better  modal  discrimination 
via  rt.  For  example,  a  1  pm  radius  device  has  ATt  =  T^HEMll)  —  rt(TE01)  =  0.0759 
for  the  etched-post  VCSEL,  compared  to  only  0.0379  for  the  oxide-apertured  VCSEL.  This 
illustrates  the  effectiveness  of  employing  a  small  radial  index  difference  to  introduce  mode 
selective  losses  and  enhance  single  mode  lasing  [25]. 

In  Figs.  13  and  14  we  plot  the  resonant  wavelengths  as  a  function  of  cavity  radius  for 
various  modes  in  etched-post  and  oxide-apertured  structures.  The  resonant  wavelength  blue- 
shifts  as  the  oxide  or  post  diameter  shrinks,  a  dramatic  departure  from  plane-wave  results. 
This  wavelength  shift  can  be  easily  explained  by  examining  the  weighted  dispersion  relations 
(21)  and  (24),  both  of  which  take  the  functional  form 

<*;>*<#>  =  ^<S>-  (75) 

Although  both  (k2)  and  (/32)  change  slightly  as  the  radius  shrinks,  {(32)  remains  very  close 
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FIG.  8:  Index  and  standing  intensity  profile  along  the  axial  direction  for  a  1.4  /j,m  radius  oxide- 
apertured  VCSEL. 

to  its  plane-wave  value.  Therefore,  as  (k2)  increases  from  its  plane-wave  value  of  zero,  u 
increases,  leading  to  the  blue  shift  illustrated  in  Figs.  13  and  14  .  This  effect  has  been 
previously  estimated  using  simpler  approximations  [26] ,  but  the  self-consistent  WIM  results 
should  be  more  accurate.  To  construct  low-threshold  microcavity  VCSELs,  the  quantum 
well  emission  peak  must  be  matched  to  the  blue-shifted  cavity  resonance  of  the  desired  lasing 
mode. 

The  resonant  wavelength  changes  more  quickly  with  radius  in  oxide-apertured  VCSELs 
as  compared  with  the  etched-post  structures.  This  occurs  because  the  oxide-apertured 
device  exhibits  a  larger  difference  in  effective  index  between  the  inner  and  outer  radial 
regions,  resulting  in  a  larger  field  confinement  and  a  correspondingly  larger  value  for  (k2) 
in  the  waveguide  core.  As  a  result,  the  oxide-apertured  structure  provides  more  spectral 
discrimination  between  the  VCSEL  resonant  modes.  For  example,  in  a  1  fim  VCSEL  AA  = 
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FIG.  9:  HEM11  mode  standing  field  intensity  for  a  1.4  pm  radius  etched-post  VCSEL.  The  field 
intensity  on  the  top  surface  is  amplified  in  order  to  illustrated  the  emitted  mode. 

A(HEMll)  —  A(TE01)  =  68.9  A  for  the  oxide-apertured  structure  as  compared  to  59.8  A 
for  the  etched-post.  This  effect  might  be  exploited  by  tailoring  the  spontaneous  emission  of 
these  devices  to  create  lower  threshold,  higher  efficiency  devices. 

Finally,  we  plot  the  components  of  the  quasi-degenerate  LPn  mode  in  Figure  15.  The 
results  show  a  small,  but  non-zero,  splitting  of  the  mode  for  aperture  radii  <~  0.85  pm, 
indicating  where  the  LP  mode  approximation  breaks  down. 

Within  the  WIM  framework,  we  approximate  the  VCSEL  cavity  modes  as  superpositions 
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FIG.  10:  TE01  mode  standing  field  intensity  for  a  1.4  pm  radius  oxide-apertured  YCSEL.  The 
field  intensity  of  the  top  surface  is  amplified  in  order  to  illustrated  the  emitted  mode. 

of  cylindrical  waves.  Each  of  these,  in  turn,  can  be  viewed  as  a  superposition  of  TE  and  TM 
plane  waves  propagating  at  an  angle 

(0)  =  arctan  (76) 

to  the  ^  axis.  As  the  cavity  radius  decreases,  the  effective  transverse  propagation  constant 
(kp)  increases  for  all  modes  and,  consequently,  the  average  angle  of  incidence  for  the  compo¬ 
nent  plane  waves  impinging  on  each  DBR  interface  increases.  The  power  reflectivity  for  TE 
waves  increases  monotonically  with  angle  until  the  total  internal  reflection  angle  is  reached, 
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FIG.  11:  Transverse  confinement  factor  for  the  first  two  modes  of  the  etched-post  and  oxide- 
apertured  VCSEL. 


FIG.  12:  Transverse  confinement  factor  for  the  fundamental,  and  a  few  sample  higher  order  modes 
for  the  oxide-apertured  VCSEL. 


while  the  power  reflectivity  for  TM  waves  decreases  monotonically  until  the  Brewster  angle 
is  reached  [27].  As  a  result,  the  TE  wave  components  encounter  more  reflective  DBR  mirrors 
as  the  cavity  radius  shrinks,  while  the  TM  wave  components  encounter  less  reflective  DBRs: 
Fig.  16  illustrates  this  behavior  for  the  oxide-apertured  structure.  Therefore  the  TE  wave 
components  require  less  threshold  material  gain,  the  TM  wave  components  more,  and  the 
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FIG.  13:  Resonant  wavelength  for  the  first  two  modes  of  the  etched-post  and  oxide-apertured 
VCSEL. 


FIG.  14:  Resonant  wavelength  for  the  fundamental,  and  a  few  sample  higher  order  modes  for  the 
oxide-apertured  VCSEL. 

threshold  gain  for  pure  TE  VCSEL  modes  decreases  with  cavity  radius,  while  the  threshold 
gain  for  TM  modes  increases. 

HEM  modes  contain  both  TE  and  TM  wave  components,  approximately  canceling  out 
the  changing  reflectivity  effects.  All  of  these  trends  are  evident  in  Figs.  17-18,  which  show 
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FIG.  15:  Resonant  wavelength  for  the  components  of  the  quasi-degenerate  LP11  mode,  illustrating 
the  point  at  which  the  degeneracy  is  broken. 


FIG.  16:  Top  DBR  reflectance  for  the  TE  and  TM  modes  of  the  oxide-apertured  VCSEL. 
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FIG.  17:  Theshold  gain  for  the  fundamental,  and  transverse  electric  and  magnetic  modes  for  the 
etched-post  YCSEL. 

the  threshold  material  gain  gth  =  AitKactive/\  for  the  fundamental  and  several  higher  order 
VCSEL  transverse  modes.  These  results  agree  with  prior  in-plane  laser  studies  that  show  a 
higher  facet  reflectance  for  TE  than  TM  modes  and  a  similar  propagation  constant  (f3)  for 
both  TE  and  TM  modes  [17,  28-30].  For  both  VCSELs  and  in-plane  lasers,  these  trends  in 
mirror  reflectivity  and  the  associated  lasing  threshold  become  more  pronounced  for  higher- 
order  TE  modes. 

Based  on  the  above  arguments,  we  expect  the  TE  mode  thresholds  to  decrease  monoton- 
ically  with  cavity  radius  and  the  TM  and  HEM  thresholds  to  increase  monotonically.  The 
TM  and  HEM  modes  generally  behave  as  expected  in  Figs.  17  and  19,  but  the  TE  mode 
thresholds  increase  abruptly  below  a  critical  cavity  radius.  This  increase  results  from  a  dif¬ 
ferent  effect — the  loss  of  mode  confinement  as  the  cavity  shrinks.  Recall  that  the  transverse 
confinement  factor  R  determines  the  strength  of  the  coupling  between  a  given  cavity  mode 
and  the  active  quantum  well,  and  that  R  decreases  monotonically  with  cavity  radius  for  all 
modes.  As  R  decreases,  the  material  gain  must  increase  accordingly  for  the  cavity  mode  to 
reach  threshold.  At  fairly  large  cavity  radii  this  effect  is  dominated  by  the  change  in  DBR 
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FIG.  18:  Theshold  gain  for  the  fundamental,  and  a  few  sample  higher  order  hybrid  modes  for  the 
oxide-apertured  VCSEL. 

reflectivity,  and  the  TE  mode  thresholds  continue  to  decrease  with  cavity  radius.  But  below 
a  critical  cavity  radius  the  modal  confinement  effect  dominates,  and  the  quantum  well  gain 
must  increase  sharply  for  the  cavity  mode  to  reach  threshold.  Eventually  the  transverse 
confinement  becomes  so  weak  that  the  cavity  cannot  support  a  radially-bound  mode  at  all. 
The  radii  at  which  our  curves  terminate  in  Figs.  17-18  indicate  this  “minimum”  cavity  size 
for  each  optical  mode.  This  size  effect  has  been  previously  estimated  by  simpler  calculations 
[26] ,  but  the  WIM  values  should  be  more  accurate.  Since  our  oxide-apertured  structure  con¬ 
fines  the  fields  better  than  the  etched-post  structure,  it  supports  bound  modes  at  smaller 
cavity  radii. 

To  summarize  the  WIMP  formalism,  we  have  presented  an  extremely  general  and  rapid 
technique  for  estimating  the  spatial  profile,  optical  confinement  factor,  resonant  frequency, 
and  threshold  material  gain  of  lasing  modes  in  cylindrically- symmetric  VCSEL  geometries. 
Variational  principles  dictate  that  this  method  will  generate  the  best  solution  over  the  space 
of  all  separable  functions.  We  also  expect,  on  variational  grounds,  that  the  “eigenvalues” 
corresponding  to  each  mode — the  resonant  frequency  and  threshold  gain — will  be  more  ac- 
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FIG.  19:  Theshold  gain  for  the  fundamental,  and  transverse  electric  and  magnetic  modes  for  the 
oxide-apertured  VCSEL. 

curate  than  the  wavefunction  itself.  Thus,  among  separable  or  single-term  approximations 
to  the  lasing  modes  this  method  should  yield  the  best  possible  results.  The  fundamental 
limitation  to  this  technique  is  the  base  assumption  of  separability.  Therefore  non-separable 
effects,  such  as  mode  divergence  and  possible  parasitic  diffraction  losses,  are  not  addressed. 

We  have  applied  the  technique  to  oxide-apertured  and  etched-post  VCSELs,  predicting  a 
blue-shift  of  the  lasing  wavelength  as  the  cavity  radius  shrinks.  It  is  essential  to  incorporate 
this  blue  shift  when  designing  the  optical  cavity,  since  the  lasing  resonance  must  line  up 
precisely  with  the  quantum  well  gain  peak  in  order  to  minimize  the  threshold  current.  Our 
calculations  suggest  that  quantum  well  emission  should  be  tuned  to  the  first  TE  mode 
cavity  resonance,  since  this  mode  exhibits  a  decreasing  threshold  gain  as  the  cavity  radius 
shrinks.  However,  below  a  critical  cavity  radius  the  threshold  gain  increases  rapidly  and 
lasing  becomes  impossible.  Our  model  can  be  used  to  investigate  a  variety  of  cavity  designs 
in  an  effort  to  minimize  this  critical  radius,  and  thereby  minimize  the  VCSEL  threshold 
current.  In  order  to  actually  calculate  this  threshold  current ,  it  will  be  necessary  to  merge 
this  model  with  quantum  well  gain  and  emission  models. 
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WIMP  Technology  Transition  Note 


The  WIMP  simulation  source  code,  along  with  bandstructure  source  code  for  quantum 
well  gain  calculations  and  associated  source  code  for  predicting  threshold  currents  as  a  func¬ 
tion  of  oxide  aperture  diameter  and  other  relevant  device  parameters,  was  transferred  to  CFD 
Research  Corporation  ( CFDRC),  of  Huntsville,  Alabama.  The  vehicle  for  this  technology 
transfer  was  a  Cooperative  Research  and  Development  Agreement  (CRADA),  formalized  by 
signatures  from  both  parties  in  October,  2000.  Under  this  agreement,  AFRL/SND  obtained 
in  return  a  licensed  set  of  CFDRCs  bundled  suite  of  simulation  tools  for  semiconductor  de¬ 
vices  (electrical,  thermal,  and  optical  modelling,  as  well  as  some  MEMS  capability).  Since 
this  time,  CFDRC  was  acquired  by  ESI  Group  on  January  28,  2004.  At  the  time  of  this 
writing,  many  issues  regarding  customer  support,  DARPA  commitments,  and  other  respon¬ 
sibilities  of  CFDRC  were  still  being  settled  as  a  result  of  this  merger.  The  ESI  Group  is 
headquartered  in  Paris,  with  U.S.  offices  Huntsville,  San  Diego,  and  Columbia,  Maryland. 


Vector  Finite  Element  Method  (VFEM) 


In  this  section,  we  present  a  quasi-exact  finite  element  method  (FEM)  model  for  analyzing 
the  optical  modes  of  microcavity  VCSELs.  This  model,  like  the  WIMP  code  describe  in  the 
previous  section,  is  also  based  on  a  variational  solution  of  the  vector  Helmholtz  equation  in 
microcavity  geometries.  The  results  of  this  model  not  only  allow  for  direct  calculation  of 
lasing  mode  parameters,  but  also  a  better  understanding  of  the  underlying  physics  associated 
with  VCSEL  oxides. 

To  begin,  we  want  to  find  the  electric  and  magnetic  field  prohles,  the  resonant 
wavelength  ( A0),  and  the  threshold  gain  (gth)  f°r  each  cavity  mode  in  azimuthally-symmetric 
VCSEL  structures.  For  this  we  must  solve  Maxwell’s  equations  subject  to  appropriate 
boundary  conditions  at  each  material  interface.  The  steady-state,  time-harmonic  electric 
held  satisfies  the  vector  Helmholtz  equation  (in  MKS  units) 


V  x 


— V  x  E 

hr 


(77) 


where  E,  H,  and  the  electric  current  J  depend  on  time  as  el,Jjt  {u  =  27tc/A);  e0)  bo •  and  k0 
are  the  free  space  permittivity,  permeability,  and  propagation  constant,  respectively.  It  can 
be  shown  [31]  that  a  weak  solution  to  (77)  may  be  obtained  by  extremizing  the  functional 
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(By  extremize,  we  mean  forcing  <51  =  0,  where  S  is  the  variational  operator.)  Here  Q  and  T 
are  the  problem  domain  and  boundary,  respectively  For  azimuthally-symmetric  structures, 
the  material  parameters  ( pr  and  er  )  are  functions  of  p  and  2  only,  and  we  may  separate  out 
the  0  dependence  in  (78)  by  assuming 

E<f,  ^  cos  (m0) ,  (79) 


Ep,  Ez^  sin  (m0) ,  (80) 

where  m  is  the  azimuthal  mode  number.  For  integer  m,  the  integrals  over  0  =  [0,  27t]  yield 
a  constant  factor  which  may  be  ignored,  effectively  reducing  the  dimension  of  the  problem 
from  three  to  two.  In  addition,  since  different  azimuthal  modes  are  orthogonal,  we  may  deal 
with  each  value  of  m  independently. 

For  lasing  mode  analysis,  we  set  the  source  current  J  to  zero,  making  (77)  a  source-free 
eigenmode  problem.  By  assuming  perfect  conducting  boundary  conditions  (n  •  E  =  0^  on 
T,  which  we  justify  later,  the  surface  integral  in  (78)  drops  out  and  we  are  left  with 


V  x  E)  ■  (V  x  E)  dv  —  kl  /  /  erE  ■  E  dv. 


Here  fl  represents  the  two-dimensional  domain  (of  the  VCSEL)  over  the  p  —  z  plane.  Our 
task  is  then  to  find  the  E  field  which  extremizes  (81).  That  is,  we  must  find  the  E  that 
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for  all  field  variations  8E.  Equation  (82)  is  a  generalized  eigenvalue  problem, 

S{E)-£T(E)  =  0. 


The  S  and  T  operators  are  defined  by 


5(,)  s  yy0  x  •  iv  x  *’■ 
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and 


(85) 


The  eigenvalue  £  is  defined  as 


T(T)  =  [[  er  T  •  6$  dv. 

J  Jvl 


F  =  k2-  — 
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All  the  desired  mode  information  may  be  found  by  solving  (82).  For  each  mode,  the  eigen¬ 
value  is  the  square  of  the  (generally  complex)  free-space  propagation  constant,  which  is 
related  to  the  modal  wavelength  (Ao)  (98)  and  the  total  optical  loss  or  threshold  gain  (gth) 
(102).  The  eigenvectors  are  simply  the  time-harmonic  (vector)  electric  fields  .  We 
approximate  the  solution  to  (82)  using  the  finite  element  method. 

The  formalism  developed  thus  far  is  general  and  could  apply  to  any  of  several  variational 


approaches  to  solving  Maxwell’s  equations.  We  now  narrow  our  attention  to  the  finite 
element  method  (FEM).  In  the  finite  element  method,  the  solution  to  (82)  is  approximated 
by  limiting  the  space  of  admissible  functions  E  to  the  linear  superposition  of  a  finite  set 
of  basis  functions  (generally  characterized  by  the  fact  that  they  are  non-zero  only  over  a 
subdomain  fle,  the  domain  of  mesh  element  e).  In  our  vector  FEM,  we  expand  the  fields 
over  a  basis  of  vector  functions  [32] , 


N 

E  =  y (87) 

2=1 

Here  N  is  the  total  number  of  basis  functions  in  the  expansion,  aq  are  unknown  coeffi¬ 
cients,  and  xpi  are  the  vector  basis  functions.  These  functions  are  second  order  node,  edge, 
and  face  element  functions,  given  in  Appendix  .  They  are  specifically  designed  to  model 
m  =  1  modes.  Substituting  (87)  into  (82)  and  exchanging  the  order  of  summation  and 
integration  yields 


xt  |  JJ  —  (\/  x  'ifti'j  ■  x  dv  —  kl  JJ  cr?/y  •  ifj  dv  =  0.  (88) 

Ensuring  (88)  holds  for  all  (same  set  of  functions  as  ipf),  ensures  it  will  hold  for  any 
linear  superposition  of  xpj ,  the  finite  basis  analogy  of  5E.  The  N  equations  represented  by 
(88)  are  exactly  (83)-(85)  taken  over  a  finite  basis,  written  conveniently  in  linear  algebra 
notation  as 

SX— £TX  =  0,  (89) 
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FIG.  20:  Sample  finite  element  mesh  for  an  oxide-apertured,  oxide  DBR  VCSEL. 

where, 

si,j  =  JJ  ““  (V  x  fi)  •  (V  x  $i)  dv »  (90) 

U,j  =  /  /  £rA-  $3  dv,  (91) 

J  Jn 

and  the  eigenvalue  (0  definition  (86)  remains  unchanged. 

We  define  our  basis  set  over  a  triangular  mesh  in  the  p  —  z  plane,  as  illustrated  in  Fig. 
20.  By  using  (randomly  shaped)  triangular  elements,  we  can  accommodate  general  VCSEL 
designs — including  for  example  tapered  oxides — and  avoid  creating  an  artificial,  mesh-driven 
predisposition  to  any  given  vector  field  component.  Due  to  the  form  of  our  basis  expansion 
and  the  use  of  absorbing  regions  (discussed  below),  S  and  T  will  be  very  large  (. N  ~  50,  000), 
sparse,  non-Hermitian  matrices.  As  a  result,  special  matrix  techniques  are  required  to  solve 
(89).  In  practice  this  is  the  most  challenging  part  of  the  finite  element  solution,  and  certainly 
the  most  time  consuming.  We  chose  to  solve  the  eigenvalue  problem  using  an  iterative 
Arnoldi  algorithm  [33,  34]  with  spectral  transformation;  this  algorithm  allows  us  to  search 
for  “mildly  complex”  (e.g.,  Re(£)  Im(£))  eigenvalues  over  a  given  range  of  the  real  axis. 


41 


By  far  the  greatest  potential  error  source  in  the  finite  element  VCSEL  analysis  is  mesh 
termination.  Due  to  the  incomplete  optical  confinement  of  both  the  oxide  apertures  and 
the  VCSEL  mirrors,  the  true  domain  of  a  VCSEL  field  solution  extends  to  infinity  in  all 
directions.  Hence,  an  artificial  mesh  termination  is  required  for  FEM  application.  A  prop¬ 
erly  designed  mesh  termination  for  this  problem  must  mimic  the  unbounded  nature  of  the 
domain,  eliminating  non-physical  reflections  at  the  mesh  edge.  In  addition,  since  the  FEM 
computational  demand — both  CPU  time  and  memory — scales  super-linearly  with  N,  we  pre¬ 
fer  to  place  our  termination  as  close  as  possible  to  the  primary  domain  of  interest,  thereby 
minimizing  the  amount  of  “wasted”  mesh  space. 

To  terminate  our  mesh,  we  insert  an  artificial  absorbing  layer  (AL)  between  the  principle 
problem  (VCSEL)  domain  (tty)  and  the  problem  boundary  (r)  [31].  This  layer  allows 
us  to  use  perfect  conductor  boundary  conditions  on  T,  as  we  assumed  earlier,  eliminating 
the  surface  integral  term  in  the  variational  form  and  dramatically  simplifying  the  FEM 
analysis.  We  define  fly  by  a  rectangular  region  in  the  p  —  z  plane,  bounded  by  the  bottom 
mirror-to-substrate  plane  (Tb)  and  the  top  mirror-to-air  plane  (Tt)  in  the  z  direction,  and 
the  transverse  lasing  mode  size  (Ts)  in  the  p  direction.  We  determine  Tg  a-posteriori ,  as 
discussed  below.  We  surround  fly  by  the  AL  as  illustrated  in  Fig.  21,  where  the  AL  domain 
is  the  union  of  the  top,  side,  and  bottom  AL  regions  (Dal  =  DT  fl  D5  fl  D#).  Due  to  the 
high  reflectance  of  the  distributed  Bragg  reflectors  (DBRs),  our  main  concern  is  absorbing 
any  radiation  incident  on  the  radial  boundary  T^.  Therefore,  we  focus  our  analysis  on  the 
radial  AL  (D5);  the  optimal  axial  AL  design  ( and  Dg)  falls  out  of  the  radial  analysis  as 
we  show  explicitly  later. 

A  basic  requirement  for  minimizing  reflection  is  that  the  impedance  (y p/e )  of  the 
absorbing  layer  in  each  region  must  match  the  radially  adjacent  VCSEL  region.  We  enforce 
this  condition  by  defining  the  absorbing  layer  material  parameters  as  1 

£rj(p)  =  a(p)er,j,  and  Prj(p)  =  (92) 

Here  srj  and  /irj-  are  the  VCSEL  material  parameters  in  axial  region  j  (defined  by  Zj-±  < 

1  Before  choosing  (92),  we  considered  a  diagonal  anisotropic  a(p)  for  a  perfectly  matched  layer  (PML)  design 
[35,  36].  However,  we  found  that  the  extra  degrees  of  freedom  provided  no  advantage  for  absorbing  general 
outward  propagating  cylindrical  waves  of  the  form  ( H ( kp )  e'lPzezrn(f)\ 
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FIG.  21:  Illustration  of  the  FEM  YCSEL  problem  domain:  Fly  is  the  YCSEL  domain,  and  f 1$,  FIt, 
and  ily  are  the  side,  top,  and  bottom  absorbing  layer  (AL)  domains,  respectively.  The  VCSEL 
and  AL  domains  are  separated  by  the  boundaries  Ts,  Ft,  and  T b,  and  the  entire  problem  domain 
is  bounded  by  the  closed  cylinder  T. 

^  <  Zj)  radially  adjacent  to  the  AL  (e.g.,  just  to  the  left  of  Tg),  and 

a{p)  =  1  —  ib{p) .  (93) 

The  absorbing  layer  performance  depends  entirely  on  the  function  b(p).  We  arrive  at 
a  suitable  function  b(p)  through  a  detailed  design  optimization.  First,  we  obtain  a  rough 
estimate  using  an  asymptotic  description  of  the  fields  and  reflections,  as  described  in  more 
detail  later.  Then  we  fine-tune  the  layer  by  minimizing  the  exact  reflection  values  as  obtained 
by  a  rigorous  transfer  matrix  calculation. 

To  illustrate  the  application  of  the  FEM  for  VCSELs,  we  analyze  several  versions  of 
a  basic  870  nm  oxide-apertured,  oxide  DBR  VCSEL  [37].  The  VCSEL  has  five  and  a  half 
periods  of  Alo^Gao^As/AROy  in  the  bottom  DBR  and  four  in  the  top  DBR.  The  cavity  is  1A 
thick  and  contains  a  single  600  A  GaAs  bulk  gain  region  centered  between  two  Alo.3Gao.7As 
barrier  layers.  A  300  A  AlAs  layer  is  included  in  the  top  barrier  to  form  an  oxide  aperture. 
Although  an  actual  oxide  aperture  formed  in  AlAs  would  have  a  square  cross-section,  we 
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treat  it  as  circular  to  maintain  the  azimuthal  symmetry.  The  entire  structure  including  the 
GaAs  substrate  is  illustrated  in  Fig.  20.  The  versions  of  this  structure  we  examined  are: 

•  1A- 1THIN  the  structure  as  described  above, 

•  1A- 1  THICK  same  structure  with  a  600  A  thick  oxide  aperture, 

•  1A-2THICK  same  structure  with  a  600  A  thick  oxide  aperture  in  both  the  top  and 
bottom  barrier  regions,  and 

•  A/2-1THICK  same  DBRs  but  with  a  A/2  cavity  and  a  600  A  thick  oxide  aperture  in 
the  top  barrier  region. 

All  materials  indices  are  assumed  to  be  real  so  that  the  only  source  of  held  loss  is  through 
absorption  in  the  ALs.  For  this  particular  structure,  the  details  of  generating  the  absorbing 
layer  function  b(p)  will  later  be  shown. 

In  Figs.  22  and  23  we  plot  |  E,,\  and  the  time  averaged  electromagnetic  energy  density, 


1  r  2  2l 

w  =  —  £q £rj  E  +  /j,q  H  ,  (94) 

for  the  sample  1A-1THIN,  pox  =  0.4  /mi  VCSEL;  these  plots  are  representative  of  the 
general  lasing  mode  profiles  for  all  the  VCSELs  we  tested.  The  fields  are  found  by  sub¬ 
stituting  the  lasing  mode  eigenvector  (X)  into  the  held  expansion  given  in  (87).  The  E,p 
prohle  (Fig.  22)  is  similar  to  the  familiar  standing  wave  prohle  obtained  via  simple  scalar 
held  techniques,  however,  we  anticipate  our  result  is  more  accurate  due  to  the  full  vector 
solution.  The  energy  density  (94)  is  found  by  estimating  H  from  the  E  held  expansion  (87) 
and  Faraday’s  law.  Although  this  is  less  familiar  than  the  E4)  prohle,  it  is  a  more  accurate 
representation  of  the  spatial  mode  energy  distribution  throughout  the  VCSEL. 

Using  the  spatial  mode  prohle,  we  estimate  the  total  and  transverse  conhnement  factors 


as 


and 
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respectively;  fly  is  the  VCSEL  volume,  and  Qpump  and  Ggain  are  the  pumped  and  total 


volume  of  the  gain  region,  respectively.  For  purpose  of  calculation,  we  estimate  Cpump  as 
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FIG.  22:  Example  plot  of  |£^|  for  the  1A-1THIN  structure  with  oxide  aperture  radius  pox  =  0.4  pm. 
To  increase  clarity,  the  figure  domain  is  smaller  than  the  calculation  domain,  and  the  background 
intensity  has  been  set  to  white. 

the  volume  of  the  active  layer  inside  the  oxide  aperture  (e.g.,  with  p  <  pox).  The  total 
confinement  factor  represents  the  percentage  of  the  mode  energy  overlapping  the  active  gain 
region  2.  We  use  this  later  in  the  Fabry-Perot  laser  equation  to  estimate  material  threshold 
gain  from  the  total  modal  loss  (102).  The  transverse  confinement  factor  has  less  quantitative 
application  but  is  a  nice  indicator  of  how  well  the  lasing  mode  is  confined  in  the  transverse 
dimension. 

We  plot  rtr  for  each  of  the  VCSEL  structures  in  Fig.  24.  The  results  for  the  1A  cavity  are 
somewhat  intuitive:  for  any  given  radius  more  and  thicker  oxides  yield  greater  confinement. 
For  the  A/2  case,  the  oxide  aperture  does  not  overlap  well  (in  the  z  direction)  with  the 
standing  wave.  Therefore,  the  transverse  confinement  is  relatively  small,  however,  the  total 


2  This  is  a  somewhat  different  definition  than  the  more  standard  definition  using  the  standing  wave  intensity 
2^ 


E 


However,  the  two  definitions  should  yield  similar  results  in  most  cases. 
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FIG.  23:  Example  plot  of  the  stored  energy  density  w  for  the  1A  -1THIN  structure  with  oxide 
aperture  radius  pox  =  0.4  /mi.  To  increase  clarity,  the  figure  domain  is  smaller  than  the  calculation 
domain,  and  the  background  intensity  has  been  set  to  white. 

confinement  factor — due  to  the  short  cavity — is  very  large,  as  shown  in  Fig.  25.  The  high 
contrast  DBRs  allow  less  field  penetration  and  therefore  a  very  high  total  confinement  factor 
for  all  four  test  cases,  compared  to  analogous  semiconductor  DBR  VCSELs. 

Due  to  the  absorbing  layers,  the  eigenvalues  of  (86)  are  complex,  and  take  the  general 
form 

&  =  ^  =  ( n  +  i  ft)2  .  (97) 

cz 

Here  r,  and  ft  are  the  real  and  complex  parts  of  the  (total)  propagation  constant  k,  of 
mode  i.  From  (97)  we  immediately  recognize  the  mode  resonance  as 

A*  =  — .  (98) 

By  definition,  the  fields  vary  harmonically  as 

exp  ( iuj{t )  =  exp  ( ikid )  =  exp  ( zr^ct)  exp  (— ftcf) ,  (99) 
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FIG.  24:  Transverse  confinement  factor  verses  oxide  aperture  radius  for  the  fundamental  lasing 


mode.  The  lines  are  cubic  spline  fits  of  the  discrete  calculation  data. 


FIG.  25:  Total  confinement  factor  verses  oxide  aperture  radius  for  the  fundamental  lasing  mode. 
The  lines  are  cubic  spline  fits  of  the  discrete  calculation  data. 

where  the  mode  frequency  is  Re(u;j)  =  r;C  and  the  total  cavity  loss  rate  (1/seconds)  is 
li  =  qiC.  These  results  describe  a  three-dimensional  leaky  cavity,  losing  energy  at  a  rate  of 
li  .  We  would  prefer  to  express  our  results  in  terms  of  the  more  familiar  Fabry-Perot  (F-P) 


47 


laser  model.  To  do  this,  we  convert  c Oi  to  a  complex  propagation  constant 

A  =  Y^gain,  (100) 

where  *  is  the  complex  conjugate,  required  such  that  a  lossy  ujt  maps  to  a  lossy  A-  The 
imaginary  part  of  A  gives  the  field  loss  rate  (for  Im(A)  >  0)  in  cm,-1.  Since  energy  is 
proportional  to  the  square  of  the  field,  the  intensity  lost  per  unit  length  (cq)  is  given  by 

=  2  Im(A)-  (101) 


In  deriving  ot  all  we  have  done  is  convert  energy  lost  per  time  to  energy  lost  per  propagation 
length.  To  apply  this  loss  rate  to  a  “textbook”  Fabry-Perot  laser,  the  cavity  length  would 
have  to  be  adjusted  in  accordance  with  the  blueshift  [10].  Finally,  we  calculate  the  threshold 
gain  (gth)  from  the  F-P  lasing  condition, 

r  fVi  =  (102) 


In  Figs.  26  and  27  we  plot  the  modal  resonance  and  threshold  gain  as  a  function  of 
oxide  aperture  radius  for  all  four  test  VCSELs.  The  resonance  results  show  the  now  familiar 
blueshift,  with  more  and  thicker  oxide  apertures  yielding  a  larger  shift.  Interpretation  of  the 
threshold  gain  curves  is  more  complicated.  We  discuss  these  results  later  in  our  discussions 
regarding  optical  losses  and  diffraction  arising  in  this  model. 

To  better  understand  the  optical  loss  sources,  we  would  like  to  divide  the  total  optical 
mode  loss  (cq)  into  a  mirror  loss  ^afmirrorA ,  Jue  to  emission  out  the  ends  of  the  VCSEL, 
and  a  diffraction  loss  ^Q,idlffractlon)  j ,  due  to  loss  out  the  VCSEL  side.  To  do  this,  we  use 
conservation  of  energy  and  the  relationship  between  a  and  the  radiated  power  and  stored 
energy  [24], 


Pv 

VpWy 


(103) 


Here  vp  is  the  field  phase  velocity, 


Re 


(e  x  H'\ 


(104) 


is  the  total  time  averaged  (real)  power  exiting  the  VCSEL  (fly)  through  Ty,  and 


Wv  =  Wy  +  Wfi1 


(105) 
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FIG.  26: 


Lasing  mode  resonance  verses  oxide  aperture  radius.  The  lines  are  cubic  spline  fits  of  the 


discrete  calculation  data. 


FIG.  27:  Threshold  (material)  gain  verses  oxide  aperture  radius.  The  lines  are  cubic  spline  fits  of 
the  discrete  calculation  data. 


is  the  total  time  averaged  stored  energy.  This  consists  of  the  electric  and  magnetic  field 
energies, 


w^  = 


SqStjQ' 


E 


dv, 


(106) 


49 


and 


w$  = 


1 


IMP 


H 


dv, 


(107) 


respectively.  Equation  (103)  is  derived  from  the  relationship  between  cavity  Q  and  loss  rate, 
adapted  to  the  F-P  laser  model.  We  find  in  practice  that  numerical  application  of  (103) 
gives  results  exactly  matching  those  obtained  with  (101),  verifying  the  general  approach.  It 
is  then  a  simple  matter  to  break  a  into  radial  and  axial  components, 


and 


Here 


(diffraction) 

~  vpWv' 

(108) 

(mirror)  PT  +  Pb 

'  vpWv  ' 

(109) 

Pv  Ps  +  Pt  +  Pb, 

(110) 

breaks  the  total  power  exiting  the  VCSEL  into  that  leaving  through  the  side  ( Ps ),  and  that 
leaving  through  the  top  and  bottom  surfaces  (Pt)  and  (Pb),  respectively.  These  powers  are 
calculated  by  restricting  the  integral  in  (104)  to  the  corresponding  surfaces.  Note  that  for 
these  calculations,  we  assume  the  VCSEL  materials  are  lossless  and  calculate  “cold  cavity” 
radiative  loss  parameters.  Although  absorptive  loss  may  be  significant  in  operating  VCSELs, 
this  calculation  should  give  a  good  estimate  of  the  diffraction  loss.  The  calculation  can  be 
easily  modified  to  fold  in  absorptive  losses  if  required. 

In  Figs.  28  and  29,  we  plot  the  radial  mode  loss  and  the  radial  percentage  of  the  total 
mode  loss,  respectively.  Comparing  the  two  figures,  we  see  that  o;(mirror)  the  majority  of 
the  change  in  threshold  with  aperture  radius  is  due  to  a  change  in  a('llffract,,on).  This  is  not 
surprising,  since  we  expect  the  change  in  mirror  reflectivity  to  be  relatively  small  and  the 
diffraction  to  increase  as  the  aperture  size  decreases. 

By  examining  the  loss  characteristics  for  each  VCSEL  design,  we  may  deduce  the  physical 
mechanisms  governing  diffraction.  Our  results  support  the  idea  that  diffraction  may  be 
viewed  as  a  coupling  between  the  bound  eigenmode  and  the  continuum  of  parasitic  modes 
[38].  These  are  radially  propagating  slab  modes  in  the  unapertured  (cladding)  region. 
Deppe  explains  in  [26]  that  for  a  cavity  bound  by  perfectly  conducting  mirrors,  the  parasitic 
mode  density  will  resemble  a  slanted  staircase  following  the  three  dimensional  density  of 
modes.  The  jumps  in  the  staircase  occur  at  each  vertical  resonance  in  the  cladding  region. 
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FIG.  28:  Radial  mode  loss  verses  oxide  aperture  radius  for  the  fundamental  lasing  mode.  The  lines 


are  cubic  spline  fits  of  the  discrete  calculation  data. 


FIG.  29:  Percentage  of  the  total  mode  loss  due  to  radial  losses  for  the  fundamental  lasing  mode. 
The  lines  are  cubic  spline  fits  of  the  discrete  calculation  data. 

Deppe’s  results  suggest  that  shorter  cavities  are  superior  since  the  eigenmode  overlaps  with 
a  smaller  density  of  parasitic  modes.  This  idea  is  illustrated  in  Fig.  30. 

For  our  high  contrast  DBRs,  we  expect  very  similar  behavior  to  this  ideal  mirror  case. 
Comparing  the  threshold  gains  of  each  structure  (Fig.  27),  the  A/2  cavity  threshold  is  con- 
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FIG.  30:  Illustration  of  the  “slanted  staircase”  parasitic  mode  density.  The  density  follows  the 
three-dimensional,  free-space  density  of  optical  modes  (dotted  line).  The  jumps  in  the  staircase 
occur  at  each  vertical  resonance  in  the  cladding  region.  The  two  discrete  modes  are  representative 
of  a  low  threshold  (A),  and  a  higher  threshold  (B)  lasing  mode. 

sistently  lower  than  all  the  1A  cases.  Moreover,  comparing  the  percentage  of  the  total  mode 
energy  lost  to  diffraction  (Fig.  29),  the  A/2  VCSEL  is  again  superior.  These  results  can’t  be 
attributed  to  transverse  confinement,  since  the  A/2  VCSEL  has  the  smallest  Ttr  of  all  four 
cases.  However,  these  results  could  be  attributed  to  the  larger  total  confinement  factor  (Fig. 
25),  or  to  a  smaller  density  of  parasitic  modes  interacting  with  the  A/2  cavity  eigenmode.  It 
turns  out  that  these  two  factors  are  closely  related:  the  longitudinal  confinement  factor — or 
effective  cavity  length  in  the  cladding  region — determines  the  location  of  the  steps  in  the 
parasitic  mode  density  (Fig.  30). 

Adopting  Deppe’s  parasitic  mode  density  interpretation,  we  can  explain  the  difference 
in  threshold  between  the  A/2  and  1A  cavity  VCSELs,  but  we  cannot  easily  explain  the 
disparity  among  the  three  1A  cases,  which  should  have  very  similar  parasitic  mode  densities. 
Furthermore,  the  parasitic  mode  density  is  solely  a  function  of  (effective)  cavity  length, 
and  therefore  cannot  address  radius-dependent  changes  in  diffraction.  To  capture  these 
radius-dependent  effects,  and  to  distinguish  the  various  1A  cavity  structures,  we  propose 
that  the  diffractive  loss  is  a  function  of  both  the  density  of  parasitic  modes  and  the  coupling 
strength  between  the  eigenmode  and  the  parasitic  modes.  Moreover,  this  coupling  strength 
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is  a  function  of  two  factors:  (1)  the  relative  alignment  of  the  eigenmode  and  parasitic  mode 
propagation  vectors,  and  (2)  the  eigenmode  penetration  into  the  cladding  region. 

The  1A-1THIN  VCSEL  exhibits  the  lowest  threshold  and  the  lowest  diffraction  loss  (Fig. 
28)  of  all  three  1A  VCSELs.  Because  of  the  weak  transverse  confinement  (Fig.  24  )  and 
the  resulting  mode  spread,  the  field  in  this  structure  is  more  planar  and  the  propagation 
vector  is  paraxial  with  respect  to  normal  to  the  planar  interfaces.  As  a  result,  the  eigenmode 
wavevector  is  nearly  orthogonal  to  the  parasitic  mode  wavevectors,  which  he  principally  in 
the  plane  of  the  layer  structure.  The  parasitic  mode  wavevectors  can  have  components  out 
of  the  plane,  but  the  larger  these  components  are  the  less  energy  they  will  carry  away  in  the 
radial  direction.  We  attribute  the  low  threshold  for  this  design  to  this  misalignment  between 
the  eigenmode  and  parasitic  mode  wavevectors.  This  hypothesis  is  further  supported  by  the 
peculiar  results  for  the  double-apertured  1A-2THICK  VCSEL. 

The  double-apertured  VCSEL  has  the  largest  transverse  confinement  factor  of  all  three 
1A  cavity  VCSELs.  Using  the  weighted  index  method  [39],  we  showed  that  this  transverse 
confinement  has  two  primary  effects: 

1.  It  confines  the  mode  energy  within  the  aperture  region. 

2.  It  causes  the  eigenmode  propagation  vector  to  tilt  away  from  normal  to  the  DBRs. 

The  first  effect  acts  to  decrease  the  parasitic  mode  (or  diffractive)  loss  by  containing  the 
mode  energy,  while  the  second  acts  to  increase  it  through  stronger  coupling  (better  propa¬ 
gation  vector  alignment)  to  the  parasitic  modes.  These  processes  compete,  and,  depending 
on  VCSEL  design  and  aperture  radius,  either  effect  may  dominate.  For  the  1A-2THICK 
VCSEL,  the  transverse  confinement  is  >  0.95  for  aperture  radii  from  1.0  to  0.6  /mi;  for  these 
radii,  the  first  effect  dominates.  However,  somewhere  between  0.6  and  0.4  jam  enough  mode 
energy  exists  outside  the  oxide  aperture  such  that  the  strong  parasitic  mode  coupling  causes 
the  optical  loss  to  rapidly  increase.  In  essence,  the  cavity  in  the  unapertured  region  (p  >  pox ) 
forms  a  waveguide  whose  source  is  the  lasing  eigenmode.  When  the  eigenmode  penetration 
into  the  waveguide  region  becomes  large  enough,  the  waveguide  appears  to  “siphon”  energy 
away  from  the  eigenmode.  This  effect  is  clearly  illustrated  in  the  plots  of  \Efi\  for  the  0.6 
and  0.4  pm  cases  in  Figs.  31  and  32,  respectively. 

We  have  performed  quasi-exact  calculations  of  the  lasing  mode  parameters  in  an  oxide 
apertured,  oxide  DBR  VCSEL.  Our  results  are  consistent  with  past  trends  from  measure- 
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FIG.  31:  \E(i\  for  the  1A-2THICK  structure  with  oxide  aperture  radius  pox  =  0.6  pm.  To  increase 
clarity,  the  figure  domain  is  smaller  than  the  calculation  domain,  and  the  background  intensity  has 
been  set  to  white. 

ment  and  calculation.  We  find  a  blueshift  in  modal  resonance  with  aperture  radius  and 
design.  Most  importantly,  we  calculate  the  total  optical  mode  loss  and  threshold  gain  for 
our  microcavity  VCSELs.  From  these  results  we  found  that  the  physical  factors  which  de¬ 
termine  the  parasitic  mode  loss  are  the  density  of  parasitic  modes  and  the  percentage  of 
the  eigenmode  coupling  to  the  parasitic  modes.  The  coupling  percentage  is  itself  a  function 
of  the  eigenmode  penetration  into  the  cladding  (p  >  pox)  region  and  the  relative  alignment 
between  the  parasitic  mode  and  eigenmode  propagation  vectors.  Roughly  speaking,  these 
two  factors  determine  the  magnitude  of  the  parasitic  mode  source  and  the  strength  of  the 
eigenmode-to-parasitic  mode  coupling. 

This  model  may  be  used  to  qualitatively  explain  past,  seemingly  conflicting,  oxide  design 
results.  In  [40],  the  authors  suggest  thinner  oxides  placed  at  standing  wave  antinodes  give  the 
lowest  thresholds.  They  consider  semiconductor  DBR  VCSELs  with  longer  effective  cavity 
lengths  and  relatively  large  parasitic  mode  densities,  approaching  the  three  dimensional 
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2.5 


FIG.  32:  \E(f>\  for  the  1A  2THICK  structure  with  oxide  aperture  radius  pox  =  0.4  /im.  Note  the 
field  leakage  into  the  waveguide  formed  by  the  oxide  apertures  and  the  cavity.  To  increase  clarity, 
the  figure  domain  is  smaller  than  the  calculation  domain,  and  the  background  intensity  has  been 
set  to  white. 

density  of  modes.  As  a  result,  even  weak  coupling  to  the  parasitic  mode  continuum  will 
yield  a  large  diffractive  loss,  implying  thinner  oxides  are  better  for  their  structures.  On  the 
other  hand,  the  authors  of  [41]  work  with  a  A/2  cavity  VCSEL  with  one  semiconductor  and 
one  dielectric  DBR.  In  their  case,  the  cavity  and  corresponding  optical  density  of  modes  is 
much  smaller,  and  the  confinement  aspect  of  their  oxide  aperture  dominates  the  increased 
coupling  strength  to  yield  a  low  threshold  design. 

Finally,  we  note  that  our  general  finite  element  method  can  also  be  extended  to  calculate 
the  heat  and  carrier  distribution  throughout  the  VCSEL  to  yield  an  even  more  complete 
model.  These  are  both  scalar  quantities  and  are  therefore  much  easier  to  calculate.  The 
main  drawback  with  the  FEM  is  the  computational  power,  time,  and  memory  required  to 
solve  the  large  sparse,  complex,  non-Hermitian  generalized  eigenvalue  problem.  However, 
with  the  availability  of  much  faster  processors  and  much  larger  memory  modules  in  the  near 
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NO 


FIG.  33:  Numbering  convention  of  the  six  node,  six  edge,  and  two  face  based  element  functions. 


future,  this  issue  may  soon  disappear. 


Vector  Basis  Functions  for  m=l  Modes 


In  this  separate  section,  we  review  the  vector  basis  functions  employed  in  our  finite 
element  model.  Based  on  the  literature,  we  choose  to  span  the  mesh  using  second  order 
node,  edge,  and  face  element  functions  [42,  43].  These  functions  have  been  shown  to  be  a 
good  compromise  between  mesh  density  and  function  order,  roughly  minimizing  the  total 
number  of  unknowns  required  to  obtain  a  given  solution  accuracy.  For  each  triangle  in  the 
mesh  there  are  14  basis  functions:  6  node  based  (^N0  —  N5  j ,  6  edge  based  [Wq  —  W^j ,  and 
two  face  based  (w$,  ,  as  illustrated  in  Fig.  33[43,  44]  .  To  define  the  14  element 

functions,  we  use  simplex  (or  barycentric)  coordinates,  defined  over  the  triangular  element 
via  the  affine  transformation  [45,  46], 


1 

1  1  1 

Co 

p 

= 

Po  Pi  P2 

Cl 

z 

Zo  Zi  z2 

1 

(111) 


Here  (po,  zo),  (pi,  z\),  and  (p2,  z2)  are  the  three  corners  of  the  triangle,  and  Co,  Ci,and  (2 
are  the  simplex  coordinates.  Based  on  the  Q,  the  element  vector  functions  are 
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N0  =  (2Co  -  1)  Co  (j>  +  p),  (112) 

Ni  =  (2Ci  -  1)  Ci  (p  +  p),  (H3) 

N2  =  (2(2  ~  1)  C2  (^  +  p)  )  (114) 

-^3  =  4^oCi  (j>  +  1  (113) 

N4  =  4(^i  (“2  (p  +  /Pj  ,  (116) 

-W5  =  4(o£2  (p  +  pj  ,  (117) 


#0  =  pC1vTC2  =  ^{&2P  +  c2i}, 

#1  =  pC2VtCi  =  ^  {ffip  +  Cli}  , 

#2  =  pC2VtCo  =  ^{&oP  +  c0£}, 

#3  =  pCoVTC2  =  ^{62p  +  C2f}, 

#4  =  pCoVTCi  =  ^  {bip  +  Ciz}  , 

14^5  =  PCiVtCo  =  {^op  +  C0z}  , 

#6  =  p4Ci  (C2Vt£o  —  CoVtC2) 

=  ~ ^ ^  {C2  (pop  +  C0z)  —  Co  (b2p  +  c2z)}  , 

14^7  =  P  4C2  (CoVrCi  -  Cl  VrCo) 

pACo, 

=  ^  {Co  (&1P  +  c\z)  —  Ci  (bop  +  C0z)}  . 

Here 


d  ^  d  A 
V  T  =  -^-p  +  -^-2:, 
op  oz 

the  bi  and  cL  are  given  by  the  inverse  affine  transformation, 


(118) 

(119) 

(120) 
(121) 
(122) 

(123) 

(124) 

(125) 


(126) 
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FIG.  34:  Illustration  of  the  six  node  based  functions  (Nq  —  N$)  for  a  typical  element;  the  functions 
are  overlayed  on  a  triangle  outline  of  the  element. 
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z 

and  A  is  twice  the  area  of  the  triangle, 

2  22 

A  =  2J  Pi  («il  -  ^i2>  =  =  2_,  Zi°i 


(127) 


(128) 


(i,  il,  and  i2  are  modulo  3).  These  functions  are  illustrated  in  Figs.  34  and  35.  These  func¬ 
tions  have  three  important  properties  that  make  them  ideal  for  electromagnetic  calculation 
in  cylindrical  coordinates. 

First,  when  filling  the  S  and  T  matrices,  functions  based  on  the  same  node  or  edge  of 
neighboring  elements  (not  necessarily  with  the  same  element  function  number  N0  —  W-)  are 
assigned  the  same  coefficient  Xj.  As  a  result,  the  node  and  edge  functions  naturally  enforce 
tangential  field  continuity  (as  illustrated  in  Fig.  36)  without  necessarily  forcing  normal  field 
continuity ,  due  to  the  two  face  functions. 
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FIG.  35:  Illustration  of  the  six  edge  based  (Wo  —  W5)  and  the  two  face  based  (Wq  —  W7)  functions 
for  a  typical  element;  the  functions  are  overlayed  on  a  triangle  outline  of  the  element.  These 
illustrations  do  not  include  the  p  weighting  present  in  (118)  -  (125). 


EB 


FIG.  36:  Illustration  of  the  natural  tangential  continuity  between  mesh  elements.  By  assigning  the 
same  coefficient  X{  to  W\  of  element  A  and  Wo  of  element  B,  the  vector  sum  of  the  two  functions 
is  tangential  to  their  common  edge. 
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Second,  our  element  functions  are  modified  from  the  more  standard  form  [44]  to  force 
field  regularity  at  p  =  0  (the  axis  condition).  The  lasing  modes  we  seek  are  analogous  to 
the  (m  =  1)  LP0i  fiber  modes  [23],  where  the  azimuthal  field  dependence  is  sin(0)/  cos(0)  . 
The  proper  axis  condition  for  the  m  =  1  modes  is  [47] 

lim  EP  =  E4>,  (129) 

o 

and 

lim  Ez  =  0.  (130) 

o 

This  condition  is  included  in  (112)  -  (125)  by  weighting  the  “standard”  [44]  edge  and  face 
functions  by  p,  and  including  the  node  functions  as  part  of  the  p  field  expansion.  Note  that 
we  may  find  m  —  0  and  m  >  1  modes  using  a  basis  similar  to  (  112)  -  (125),  altered  to 
accommodate  the  proper  axis  condition  for  these  modes. 

Third  and  least  obvious,  these  elements  properly  model  the  null  space  of  the  curl  operator, 
which  has  been  shown  to  eliminate  the  occurrence  of  spurious  solutions  (modes)  [48,  49]. 

Absorbing  Layer  Design 

In  this  section,  we  give  the  details  of  an  absorbing  layer  design  pioneered  for  this  effort 
.  We  use  a  two  step  process  to  optimize  b(p),  the  absorbing  layer  loss  function:  First, 
we  obtain  a  rough  estimate  using  an  asymptotic  description  of  the  fields  and  reflections. 
Second,  we  fine-tune  the  layer  by  minimizing  the  exact  reflection  values  as  obtained  by  a 
rigorous  transfer  matrix  calculation.  In  both  steps  we  model  the  radial  AL  (fig)  as  a  set  of 
discrete  cylinders,  as  shown  in  Fig.  37,  approximating  a(p)  as  an  =  a((pn  +  pn-i)  /2),  where 
AL  cylinder  n  is  defined  by  [/>,,  -  i  ,  pn]  ■  To  simplify  further,  we  ignore  the  2  dependence  of 
the  material  parameters  and  perform  our  analysis  using  a  single  set  of  (er,  /ir),  the  rough 
mean  material  values  for  the  VCSEL.  This  should  be  sufficient,  since  we  are  not  concerned 
about  interfacial  reflections  between  the  various  axial  layers,  and  desire  only  to  suppress  non¬ 
physical  radial  reflections.  We  choose  b(p)  to  minimize  the  reflection  of  cylindrical  waves 
incident  on  the  radial  boundaries. 

These  waves  are  constructed  from  the  z  and  <p  field  components  given  in  each  region 
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FIG.  37:  Illustration  of  the  geometry  used  to  design  the  radial  AL.  is  broken  into  discrete 
cylinders,  each  with  a  constant  AL  parameter  an  =  a((pn  +  pn- 1)/2). 
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(131) 

(132) 

(133) 


(134) 


where  ao  =  1  (bo  —  0)  in  flv  and  we  have  assumed  non-magnetic  materials  ( p,r  =  1 ).  The 
coefficients  TJJo1  and  £To1  are  the  magnitude  for  the  outward  and  inward  propagating  waves, 
and  pol  =  TE  or  TM  labels  the  (uncoupled)  polarizations  of  the  m  =  0  modes.  It  has  been 
shown  [35]  that  absorbing  layers  that  perform  well  for  these  m  =  0  modes  will  also  work 
well  for  the  m  =  1  modes  that  we  are  interested  in. 

We  can  describe  cylindrical  waves  incident  at  an  angle  9  — measured  from  the  normal  to 
Ts  in  the  p  —  z  plane  (hps) — to  the  radial  interface  by  writing  the  radial  ( kn )  and  axial  (f3n) 
propagation  constants  as 


kn  =  k0r/an  cos  (9) , 
(3n  =  kopan  sin  (9) . 
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(135) 

(136) 


This  choice  ensures  that  kn  and  3n  satisfy  the  dispersion  relation 


uj2n0£o£ra2n  =  kl£ra2n  =  k2n  +  fa, 


in  each  region,  where  the  (generally  complex)  index  of  refraction  is 


0  =  ^£r 


(137) 


(138) 


For  the  initial  AL  design  analysis  we  employ  the  theory  of  small  reflections  [50] ,  approx¬ 
imating  the  total  Fresnel  reflection  for  a  cylindrical  wave  incident  on  the  AL  at  p0  as 


ro|  +  \ri\  h+\  r2\  *1,2  +  •  •  •  +  rNcyl 


(139) 


Here  Ncyi  is  the  total  number  of  cylinders  in  the  AL,  rn  is  the  Fresnel  reflection  coefficient 
at  pn,  and  gives  the  attenuation  due  to  the  propagation  from  p0  to  pn  and  back. 

We  use  |  rn  |  rather  than  rn  in  (139)  to  minimize  the  interference  effects,  since  we  are  after 
a  broadband  AL  design  optimized  for  all  angles  of  incidence.  We  estimate  rn  using  the 
asymptotic  form  for  cylindrical  waves  in  each  region  n.  Due  to  the  form  of  (92)  and  (131) 
(134),  the  magnitude  of  the  TE  and  TM  Fresnel  reflection  coefficients  are  the  same,  hence  we 
need  only  perform  the  analysis  once.  Taking  the  asymptotic  form  for  the  Hankel  functions, 


hm  H( j1)/(2)  (M  -  e-=, 
P^°°  VknP 


(140) 


(+  applies  to  (1)  and  —  applies  to  (2))  and  enforcing  tangential  held  continuity,  we  find 


kn-\- 1  ®n+ 1  kn 

O'nkn+l  T  &n+\kn 

COS  (0n+l)  ~  COS  (Qn) 
COS  (9n+ 1)  +  COS  (0n) 


(141) 


for  n  =  0  to  Ncyi  —  1.  6n  and  0n+\  are  the  (complex)  propagation  angles  with  respect  to 
hrs  in  regions  n  and  n  +  1,  respectively.  We  set  |  rN,:fl,  |  =  1  to  enforce  the  perfect  conductor 
boundary  condition  at  the  mesh  edge  (T).  Keeping  only  the  first  (lowest  order)  term  in  the 
derivatives  of  (140),  we  approximate 


tl,2 f if 2  •  •  •  tni 


(142) 


where 


f  ==  f:)  i2knln 
Ln  —  c 


(143) 
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and  kn  and  ln  =  pn  —  pn_i  are  the  radial  propagation  constant  and  cylinder  thickness  in 
region  n.  Substituting  (141)  -  (143)  into  (139)  and  taking  the  magnitude  yields, 


r 


NCyl  Tl 

X>n 


n= 0  n— 0 


(144) 


which  we  may  solve  for  a  given  a(p)  to  find  |  r  ( 0t )  | ,  the  magnitude  of  the  total  Fresnel 
reflection  for  a  general  (e.g.,  includes  axial  propagation  exp (i/3z))  cylindrical  wave  incident 
on  the  AL  at  angle  +  Note  that  (144)  is  equivalent  to  the  result  obtained  via  a  plane  wave, 
planar  interface  analysis,  and  may  therefore  be  directly  applied  to  the  axial  AL  design  (f+ 
and  f+),  in  addition  to  the  radial  design  (f+). 

To  finalize  our  radial  AL  design  we  use  a  2  x  2  transfer  matrix  solution  for  the  TE  or  TM 
fields.  Enforcing  continuity  of  (131)  and  (134),  or  (132)  and  (133)  at  each  radial  interface 
pn,  we  have 


L 


n 


K 


=  R 
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(145) 


where 


and 
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fff)(t+ ip„)  iffiw.) 


(146) 


(147) 


We  have  dropped  the  pol=TE/TM  superscript  since  (145)  holds  for  both  polarizations.  At 
the  mesh  boundary  (r),  we  enforce  E4)  =  0  which  relates  the  coefficients  in  the  outermost 
AL  cylinder  by 


ANCyi  ~  ~BNcyl 


H?{kNcylPNcyl) 


Hf\kNcvlPNcvl) 


(148) 


For  convenience,  we  choose  ANcyl  =  - H^\kNcylpNcyl )  and  BNcyl  =  h[2) (kNcylpNcyl).  Begin¬ 
ning  in  the  outer  AL  cylinder  and  working  inward  via  repeatedly  applying  (145),  we  calculate 
A0  and  B0.  The  total  Fresnel  reflection  coefficient  for  an  outward  propagating  wave  incident 
on  the  AL  is  given  by 


r 


Bo 

+ 


(149) 
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To  design  our  AL,  we  ran  the  analysis  outlined  previously  with  VCSEL  material  param¬ 
eters  sr  =  9  and  jir  =  1.  We  assumed  a  polynomial  form  for  the  AL  loss  function, 

order 

Kx)  =  ~hiX% ■  (15°) 

2=1 

Here  order  refers  to  the  polynomial  order,  bt  are  the  polynomial  coefficients,  and 

P  ~  Prs 

y  =  - — 

PT  ~  Pv  s 

is  the  normalized  radial  coordinate  ranging  from  0  at  Ts,  the  VCSEL  AL  interface,  to  1  at 
T,  the  problem/mesh  boundary.  Through  a  numerical  optimization,  we  found  that  a  second 
order  polynomial  with  coefficients, 

b\  =  0.1178,  (152) 

b2  =  0.7433,  (153) 

worked  best  for  b( y).  In  addition,  we  found  an  AL  thickness  of  1.5A  (at  A  =  870  nm) 
was  superior,  with  little  change  in  the  AL  properties  for  greater  thicknesses.  The  reflection 
coefficient  (calculated  with  the  transfer  matrix  approach)  verses  angle  is  given  in  Fig.  38. 
The  key  point  to  recognize  from  Fig.  38  is  the  intensity  reflection,  R,jb  =  201og10(|  r  |)  < 
—50  dB  for  0,  <  40°,  indicating  a  good  “wide-angle”  absorber  design. 

To  provide  a  more  practical  test  of  our  AL  design,  we  analyzed  the  1A  -1THIN  VCSEL 
(previously  described),  tracking  the  lasing  mode  (defined  as  the  lowest  loss  mode,  e.g.,  the 
mode  with  the  smallest  Im(£))  eigenvalue  as  a  function  of  prs  —  pox ,  which  is  the  separation 
between  the  VCSEL-AL  boundary  and  the  oxide  aperture.  We  tested  two  different  oxide 
aperture  radii,  0.5  /mi  and  1.0  /im,  using  our  1.5A  (=  0.435  /im),  second  order  polynomial  AL 
design.  The  resulting  real  and  imaginary  parts  of  the  eigenvalues  converge  for  prs  —  pox  A 
0.9  /im  (corresponding  to  ~  3A),  denoting  the  minimum  allowable  separation  between  the 
oxide  aperture  and  Ts.  It  is  interesting  to  note  that  for  the  1.0  pm  oxide  aperture,  the 
variation  in  the  eigenvalue  from  pox  =  1.2  //in  to  2.0  /an  corresponds  to  a  variation  in  the 
resonant  wavelength  of  less  than  two  angstroms,  implying  that  the  AL  works  well  even  when 
placed  very  close  to  the  aperture.  No  analogous  test  was  performed  on  the  axial  AL  design; 
we  assume  that,  based  on  the  good  results  of  the  radial  AL  design  and  the  high  reflectance 
of  the  oxide  DBRs,  the  optimized  radial  AL  design  should  work  for  the  axial  AL  as  well. 
To  perform  the  lasing  mode  analysis,  we  used  pis  =  2.0  pm  and  varied  the  oxide  aperture 
radius  ( pox )  from  0.4  p,m  to  1.0  pm. 
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FIG.  38:  Calculated  intensity  reflection  R  =  |r|2  of  the  radial  absorbing  layer  (AL)  using  the  rig¬ 
orous  transfer  matrix  approach.  The  calculation  is  based  on  a  1.5A  thick  AL  at  870  nm,  positioned 
with  pi-s  =  2.0  //,rn.  Thirty  layers  were  used  to  discretize  the  AL  loss  function.  Inset  is  the  intensity 
reflection  in  decibels  R(ib  =  201og10(|r|). 

Comparison  of  WIMP  and  VFEM  to  Published  Results 

In  this  section,  we  present  detailed  comparison  and  results  when  both  the  weighted  index 
method  and  the  vector  finite  element  method  simulations  are  applied  to  two  structures  found 
in  the  literature.  The  goal  here  is  not  only  to  calibrate  and/or  validate  each  model  to  known 
results,  but  also  to  examine  the  individual  strengths  and  shortcomings  of  each  model.  As 
a  quasi-exact  approach,  we  expect  the  VFEM  to  give  superior  results  for  all  investigated 
parameters  (resonant  wavelength,  transverse  confinement  factor,  diffraction  losses,  etc. . .). 
Nevertheless,  the  computation  time  and  problem  setup  may  make  this  method  unattractive 
for  a  number  of  applications,  and  therefore  this  comparison  of  the  two  methods  is  warranted. 

We  computed  the  resonant  wavelength,  threshold  gain,  and  transverse  confinement  fac¬ 
tor  as  a  function  of  oxide  aperture  radius  for  four  variations  of  VCSEL  structures  grown, 
fabricated,  and  tested  by  Professor  Y.-H.  Lee’s  research  group  at  the  Korea  Advanced  In¬ 
stitute  of  Science  and  Technology  (KAIST)[51,  52],  It  can  be  seen  from  Fig.  39  that  the 
eigenmode  resonances  match  very  well  for  each  structure.  The  WIMP  results  deviate  by  a 
uniform  translation  in  wavelength  by  ~  2 nm.  This  shift  may  be  due  to  the  difference  in 
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FIG.  39:  Comparison  of  the  resonant  wavelength  calculated  using  the  WIMP  and  VFEM  methods 
for  the  four  variations  of  the  KAIST  YCSEL  described  in  Refs. [51,  52], 

how  the  VCSELs  are  represented  in  this  versus  the  VFEM.  Such  small  discrepancies  can 
easily  be  attributed  in  how  the  material  layer  thicknesses  or  indices  of  refraction  are  han¬ 
dled  in  these  codes.  The  threshold  gain  (Fig. 40)  shows  much  larger  discrepancies,  however. 
Indeed,  the  WIMP  values  for  gth  deviate  by  roughly  2x  the  values  of  the  reference  VFEM 
results  for  all  three  1A  cases.  The  A/2  thick  cavity  results  are  even  worse,  with  the  wrong 
trend  even  occurring:  the  WIMP  has  the  A/2  case  having  the  highest  value  of  gth,  while 
the  VFEM  demonstrates  the  experimentally  validated  case  of  such  cavities  having  lower 
threshold.  This  lack  of  agreement  could  be  partially  explained  by  the  fact  that  this  VCSEL 
structure  stretches  the  limits  of  the  assumption  of  separability  upon  which  the  WIMP  for¬ 
malism  is  based.  More  realistically,  however,  is  a  fundamental  problem  in  how  the  diffraction 
losses  are  calculated  in  WIMP.  Finally,  the  transverse  confinement  factor,  rtr,  is  displayed 
for  these  structures  in  Fig. 41.  From  this  figure,  it  is  evident  that  both  models  agree  very 
well.  In  total,  then,  the  WIMP  is  seen  to  give  reasonable  results  (and  certainly  good  quali¬ 
tative  trends !)  for  predicting  A  and  F,  but  parasitic  mode  coupling  is  still  an  issue  of  further 
research. 
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FIG.  40:  Comparison  of  threshold  gain  calculated  using  the  WIMP  and  VFEM  methods  for  the 
four  variations  of  the  KAIST  VCSEL  described  in  Refs.  [51,  52]. 


FIG.  41:  Comparison  of  the  transverse  confinement  factor  calculated  using  the  WIMP  and  VFEM 
methods  for  the  four  variations  of  the  KAIST  VCSEL  described  in  Refs.  [51,  52]. 
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FIG.  42:  Sample  VFEM  meshing  of  the  “blunt”  aperture  region  in  an  oxide-apertured  VCSEL. 


FIG.  43:  Sample  VFEM  meshing  of  the  tapered  aperture  region  in  a  theoretical  oxide-apertured 
VCSEL. 

VFEM  Applied  to  Tapered  Oxide  Apertures 

As  shown  in  the  previous  section,  the  VFEM  has  proven  to  give  superior  results  for 
diffraction  losses  occurring  in  a  variety  of  structures  when  compared  to  the  WIMP  method. 
With  its  basis  in  the  finite  element  method,  this  VFEM  code  is  also  able  to  handle  structures 
that  WIMP  simply  cannot ,  namely  devices  having  more  structure  in  the  radial  geometry. 
One  example  of  such  a  device  is  a  tapered-oxide-aperture  VCSEL.  A  VFEM  meshing  of  a 
traditional  “blunt”  oxide  aperture  structure  is  shown  in  Fig. 42,  while  that  for  a  tapered 
structure  is  shown  in  Fig. 43.  Although  simulation  of  a  structure  like  this  might  seems 
contrived,  nothing  is  further  from  the  truth.  Indeed,  it  is  possible  to  force  tapering  of  such 
structures  simply  by  varying  the  mole  fraction  of  aluminum  in  the  aperture’s  AlxGai_xAs 
layers  from,  say,  x  =  0.95  to  x  =  1.00,  which  can  be  done  either  through  mass  flow  controller 
operations  in  the  case  of  MOCVD  growth,  or  by  digital  alloying  techniques [53]. 

The  structures  under  investigation  are  a  A/2-cavity  VCSEL  based  on  a  design  from  the 
University  of  Texas  (henceforth  denoted  “UT”)[54],  and  a  4  —  A-cavity  VCSEL  design  gen¬ 
erated  at  the  University  of  Southern  California  (denoted  hereon  as  “USC”)[55].  The  UT 
structure  has  a  single  In^Gai.^As  QW  gain  source,  and  employs  AlAs  -  GaAs  bottom  DBR 
mirrors  (26  pairs)  as  well  as  post-growth-deposited  dielectric  mirrors  (6  pairs  of  ZnSe  -  MgF). 
A  schematic  representation  of  the  layer  structure  is  shown  in  Fig. 44,  and  a  plot  of  the  the¬ 
oretical  index  of  refraction  and  standing  wave  intensity  profile  is  demonstrated  in  Fig. 45. 
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FIG.  44:  Schematic  representation  of  the  “UT”  short-cavity  VCSEL  under  investigation. 

Oxide  apertures,  nominally  396  A  thick  on  either  side  of  the  QW,  occur  at  positions  midway 
between  the  cavity  standing  wave  nodes  and  anti-nodes.  In  contrast,  the  USC  “long-cavity” 
structure  is  comprised  of  all-semiconductor  mirror  layers  surrounding  the  4-A-thick  optical 
cavity.  The  structure  has  a  single  300  A  oxide  aperture  centered  on  the  standing  wave  peak 
(anti- node),  above  the  QW  active  region.  Fig. 46  demonstrates  the  schematic  representation 
and  Fig.  47  standing  wave  fields  for  this  structure,  respectively. 

For  each  of  these  structures,  we  computed  the  transverse  confinement  factor  (Ttr),  the 
diffraction  rate,  the  field  distribution  ( E '<Q,  and  the  threshold  material  gain  (gth)  with  respect 
to  oxide  aperture  taper  length  (see  Fig.  43).  For  each  case,  this  variation  was  computed 
for  three  different  inner  oxide-aperture  radii:  1,  2,  and  3  /im.  For  each  of  these  radii,  the 
taper  length  was  then  varied  from  0  pm  (“blunt”  aperture),  to  3  pm,  in  1-pm  increments. 
Additionally,  the  oxide  aperture  thickness  was  kept  constant,  and  only  the  taper  length  was 
varied,  synonymous  with  “longer  apertures  are  sharper.” 

Results  for  calculation  of  the  transverse  confinement  factor  for  both  the  USC  and  UT 
designs  as  a  function  of  aperture  taper  length  are  shown  in  Figs.  48  and  49.  As  the  taper 
length  is  increased,  the  volume  of  low  index  material  in  the  vicinity  of  the  VCSEL  core 
decreases.  This  results  in  a  reduction  of  the  effective  index[10]  contrast  between  core  and 
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FIG.  45:  Theoretical  plot  of  the  “UT”  short-cavity  VCSEL  standing  wave  fields  and  refractive 
index  profile. 


FIG.  46:  Schematic  representation  of  the  “USC”  long-cavity  VCSEL  under  investigation. 
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FIG.  47:  Theoretical  plot  of  the  “USC”  long-cavity  YCSEL  standing  wave  fields  and  refractive 
index  profile. 


FIG.  48:  Plot  of  rtr  as  a  function  of  oxide  aperture  taper  radius  for  the  UT  structure.  Note:  This 
calculation  assumes  uniform  gain  over  the  aperture  cross-section. 

cladding  regions.  Consequently,  the  mode  spreads  out  spatially,  and  Ttr  diminishes.  In 
lesson,  the  smaller  the  oxide  aperture’s  inner  radius,  the  more  pronounced  are  the  effects  of 
giving  the  aperture  a  tapered  shape. 


71 


FIG.  49:  Plot  of  F*r  as  a  function  of  oxide  aperture  taper  radius  for  the  USC  structure.  Note: 
This  calculation  assumes  uniform  gain  over  the  aperture  cross-section. 

If  we  define  the  “diffraction  rate”  as  the  percentage  of  optical  mode  energy  that  is  scat¬ 
tered  or  emitted  radially  out  the  sie  of  the  VCSEL,  than  we  can  compute  the  relative  merits 
of  tapering  an  oxide  aperture  toward  unwanted  contributions  in  effective  diffractive  losses. 
The  diffraction  rate  as  defined  in  this  manner  is  plotted  in  Fig.  50  for  the  UT  structure,  and 
in  Fig.  51  for  the  USC  design.  It  is  seen  that  tapering  the  aperture  significantly  reduces 
diffraction  losses  in  both  designs,  and  for  all  cases  (all  values  of  inner  oxide  aperture  radius). 
Based  on  the  results  for  these  two  disparate  structures,  reasons  for  this  can  be  inferred  in 
the  following  ways.  As  mentioned  previously,  more  tapering  leads  to  great  spatial  spreading 
of  the  mode.  A  resultant  flatter  wavefront  naturally  suffers  less  diffraction.  As  discussed  by 
Hegblom  and  co-workers  [11,  56],  tapering  of  an  aperture  serves  to  give  it  a  lens- like  quality, 
which  in  turn  re-focuses  the  mode  and  counters  diffraction.  Also  per  the  Hegblom  model, 
the  symmetric  taper  in  this  structure  produces  a  linear  phase  shift  as  a  function  of  radial 
position,  giving  the  element  a  lens  with  a  linear  profile.  As  a  result,  increasing  the  taper 
length  beyond  1  pm  is  seen  to  further  reduce  diffraction  in  the  USC  structure,  but  it  has  no 
appreciable  further  effect  in  the  UT  device. 

As  the  USC  cavity  length  is  nearly  double  that  of  the  UT  design,  the  longer  focal  length 
lensing  from  longer  tapering  better  matches  the  USC  cavity,  significantly  reducing  diffraction 
losses  in  this  structure  compared  to  the  UT  design.  The  mode  of  the  UT  VCSEL  is  best 
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FIG.  50:  Plot  of  “diffraction  rate”  (see  main  text)  as  a  function  of  oxide  aperture  taper  radius  for 
the  UT  structure.  The  inset  shows  a  blow  up  for  longer  taper  lengths,  demonstrating  how  slight 
the  change  in  diffraction  rate  becomes  for  increased  taper  lengths  Note:  Each  curve  corresponds 
to  a  different  value  of  oxide  aperture  inner  radius. 


FIG.  51:  Plot  of  “diffraction  rate”  (see  main  text)  as  a  function  of  oxide  aperture  taper  radius 
for  the  USC  structure.  Note:  Each  curve  corresponds  to  a  different  value  of  oxide  aperture  inner 
radius. 

suited  for  the  1  gm  tapers,  and  further  reductions  are  slight  if  at  all. 

Perhaps  the  most  dramatic  visualization  of  the  influence  of  tapered  versus  blunt  oxide 
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FIG.  52:  Intensity  plot  of  the  field  distribution  for  the  UT  structure.  On  the  left,  the  device 
as  simulated  with  no  tapering  of  the  aperture.  On  the  right,  with  a  1  pm  taper.  Note  the  radial 
“fingers”  showing  the  effects  of  diffraction  from  the  blunt  aperture. 

apertures  can  be  seen  when  we  examine  the  electric  fields  in  such  devices.  In  Figs.  52-53 
are  plotted  the  field  distributions  of  the  optical  field  for  the  UT  and  USC  structures, 
respectively.  For  the  UT  structure,  one  sees  significant  scattering  of  the  field  from  the  oxide 
aperture,  giving  rise  to  significant  energy  losses  when  no  aperture  is  included.  Similarly,  the 
“blunt”  aperture  case  of  the  USC  structure  is  also  evident.  Indeed,  with  a  little  imagination, 
one  can  almost  visualize  the  reflection  of  this  diffracted  energy  off  of  the  top  boundary  surface 
(at  10  /im  on  the  vertical  axis)  back  into  the  VCSEL  as  a  potential  guided  mode. 

Finally,  we  show  the  threshold  material  gain  (gth)  required  to  exactly  offset  optical  losses 
in  Fig.  54  (UT  structure)  and  Fig.  55  (USC  structure).  This  includes  not  only  the  effects 
of  gain-to-mode  spatial  overlap,  but  also  takes  into  account  transverse  mode  confinement 
effects  and  diffraction  losses,  as  shown  above.  For  these  results,  we  have  assumed  a  uniform 
gain  distribution  over  the  cross-section  of  the  aperture.  This  is  not  unreasonable  in  that 
one  of  the  early  motivating  factors  for  the  use  of  such  aperture  was  in  current  confinement 
or  guiding  in  transistor  type  devices.  In  a  typical  design  mindset,  for  maximizing  mode-to- 
gain  interactions,  the  optical  field  should  be  sharply  peaked  at  the  cavity  center.  However, 
as  the  VFEM  has  demonstrated,  the  more  concentrated  this  mode  is,  the  more  it  can 
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FIG.  53:  Intensity  plot  of  the  field  distribution  for  the  USC  structure.  On  the  left,  the  device 
as  simulated  with  no  tapering  of  the  aperture.  In  the  middle,  a  1  /im  taper.  On  the  right,  a  3 
gm  taper.  Again,  the  radial  “fingers”  indicate  diffraction  from  the  blunt  aperture,  with  possible 
reflection  back  into  the  structure  at  the  semiconductor /air  interface  at  10  gm  on  the  vertical  axis. 

potentially  diffract,  and  hence  there  exists  and  underlying  overlap  versus  diffraction  tradeoff 
consideration  in  optimizing  a  given  design  concept.  Indeed,  it  is  evident  in  the  calculations 
for  gth  for  the  UT  structure  that  increasing  the  oxide  taper  eventually  leads  to  increased 
gth ,  which  is  explained  as  an  eventually  loss  of  optical  confinement  losing  out  to  diffractive 
losses  in  such  a  design. 

Generalized  Interface  Approach  for  Exact  Design  in  Graded  Interface  Distributed 
Bragg  Reflectors 

In  this  subsection,  we  introduce  the  concept  of  a  generalized  interface,  and  demonstrate 
how  such  a  construct  allows  us  to  design  high  reflectivity  distributed  Bragg  reflection  mirrors 
(DBRs)  by  exactly  phase  matching  the  reflectance  in  layers  that  include  arbitrary  grades. 
The  resulting  generalized  phase  matching  constraints  depend  upon  the  phase  angles  of  the 
reflection  and  transmission  coefficients  for  these  grade  layers,  and  we  show  the  method  for 
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FIG.  54:  Threshold  material  gain,  <7^,  for  the  UT  structure.  Each  curve  is  for  a  different  oxide 
aperture  inner  radius,  and  we  have  assumed  uniform  material  gain  distribution  over  the  oxide 
aperture. 


FIG.  55:  Threshold  material  gain,  for  the  USC  structure.  Each  curve  is  for  a  different  oxide 
aperture  inner  radius,  and  we  have  assumed  uniform  material  gain  distribution  over  the  oxide 
aperture. 

calculating  these  phases  to  arbitrary  precision  either  analytically  (for  a  certain  set  of  grade 
profiles)  or  numerically  (for  arbitrary  profiles). 

The  importance  of  this  effort  follows  from  the  fact  that  the  use  of  graded  layers  instead 
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FIG.  56:  a)  Diagram  of  a  periodic  mirror  structure  with  the  various  phase  terms  illustrated,  b) 
Coordinate  system  and  basis  functions  for  the  electric  field  for  a  section  of  the  mirror. 

of  abrupt  interfaces  has  been  demonstrated  to  significantly  reduce  the  voltage  required  to 
drive  current  through  DBR  mirrors.  In  the  question  for  low-threshold  lasers,  then,  such 
graded  layer  mirrors  are  rapidly  becoming  the  norm,  and  a  correct  understanding  of  the 
optical  principles  necessary  for  microcavity  resonance  aligned  to  the  active  region’s  material 
gain  peak  is  required  to  reach  these  ultimate  low-threshold  limits. 

For  a  multilayer  mirror  to  have  high  reflectivity,  the  individual  mirror  layers  must  be 
designed  in  such  a  manner  as  to  insure  that  the  reflections  from  every  layer  constructively 
interfere.  This  is  most  easily  realized  by  ensuring  that  the  phase  change  accrued  by  the  field 
propagating  in  a  round-trip  from  one  graded  layer  to  the  next  is  equal  to  a  multiple  of  2tt. 
Figure  56  illustrates  a  periodic  mirror  structure,  and  shows  the  phase  terms  associated  with 
the  various  layers. 

The  phase  of  the  reflection  off  the  first  graded  layer  is  denoted  by  f3\ .  The  phase  matching 
condition  for  the  second  graded  layer  is  then: 
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(3\  +  27 vk  —  cx i  +  <Ph  +  /?2  “I-  &H  +  oq  k  —  1,2,...  (154) 

Here  /32  is  the  phase  for  reflection  off  the  second  graded  layer,  ay  is  the  phase  for  trans¬ 
mission  through  the  first  graded  layer,  and  <pn  is  the  phase  for  transmission  through  the 
constant  high  index  layer.  The  phase  for  constant  high  index  (n# )  and  low  index  (rq)  layers 
of  thickness  d#  and  djr,  is  simply  given  by: 


27T  Tltf  da 


4>l  = 


27T  nL  dL 


A  A 

The  phase  matching  condition  for  the  next  graded  layer  is: 


(155) 


f3 2  +  27 tI  —  CK2  +  (pL  +  /?  1  +  O/.  +0(2  l  —  1,2...  (156) 

The  transmitted  and  reflected  phase  terms  a  and  (3  are  determined  from  the  graded  layer 
designs,  and  will  depend  in  general  on  the  index  profile  and  thickness  of  the  grades.  It  is 
possible  by  recursion  to  show  that  equations  (154)  and  (156)  ensure  phase  matching  for 
every  reflection,  and  that  they  are  sufficient  to  determine  both  cm  and  cm,  and  thus  the 
thickness,  of  both  the  high  index  and  low  index  layers  for  a  given  graded  layer  design.  The 
corresponding  (minimum)  thicknesses  are: 


dn  — 


A  (2tv  +  (3\  —  @2  —  2aq) 

47T  Tin 
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where  A  is  the  operating  wavelength.  The  designs  of  the  graded  layers  are  typically 
determined  by  the  desired  electronic  transport  through  the  mirror,  with  the  constant  index 
layers  providing  the  proper  optical  thicknesses. 

In  addition  to  designing  the  layers  inside  the  mirror  stack,  it  is  also  necessary  to  have  a 
proper  thickness  for  the  layers  that  terminate  the  mirror.  In  the  case  of  a  VCSEL,  the  mirror 
would  be  terminated  on  one  side  by  the  cavity  layer,  and  on  the  other  side  by  either  air  or  the 
substrate,  depending  on  whether  the  mirror  is  the  top  mirror  or  the  bottom  mirror.  For  a 
bottom  mirror,  the  form  of  the  phase  matching  equation  remains  the  same  as  equation  (154), 
but  the  value  of  the  reflected  phase  term  (3\  may  or  may  not  remain  the  same  depending 
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on  whether  or  not  the  high  index  mirror  layers  are  composed  of  the  same  material  as  the 
substrate.  For  a  top  mirror,  the  value  of  (3\  is  zero  (since  no  grade  is  necessary). 

For  the  cavity,  the  condition  for  phase  matching  is  that  the  round-trip  phase  of  a  wave 
resonating  in  the  cavity  is  an  integral  multiple  of  2n.  Thus,  the  phase  matching  condition 
is: 


2(t)HC  +  f3i  +  @3  =  2irm,  m  =  1,2,...  (159) 

As  mentioned  previously,  the  values  of  the  phase  terms  a  and  (3  are  determined  by  the 
particular  choice  of  the  graded  layer  index  profile  and  thickness.  A  method  for  determining 
these  values  for  linear  index  (composition)  grading  will  be  given  in  the  next  section. 

The  problem  of  graded  DBR  design  has  now  been  reduced  to  a  calculation  of  the  trans¬ 
mission  and  reflection  phase  angles  a  and  (3  for  the  graded  layers.  We  obtain  these  angles 
from  the  transmission  and  reflection  coefficients  for  propagation  in  the  three  region  case 
illustrated  in  Fig.  56(b):  the  refractive  index  is  constant  in  regions  1  and  3  with  values  n\ 
and  n3,  respectively  while  region  2  is  the  grade.  The  refractive  index  varies  approximately 
linearly  with  composition  in  the  AlAs/GaAs[57]  material  system.  If  the  desired  composition 
grade  is  linear,  then  the  index  of  refraction  of  the  grade  will  have  the  form  n  (z)  -  ;/|  i  mz 
with  z  —  0  at  the  interface  between  region  1  and  region  2  and  m  =  (ri2  —  «i)  / L.  We 
limit  ourselves  to  normal  incidence  (i.e.  propagation  along  the  z  direction).  Without  loss 
of  generality,  we  can  assume  that  the  electric  field  vector  in  regions  1  and  3  consists  of  only 
a  y  component.  In  the  constant-index  regions,  the  field  solutions  are  simply  plane  waves 
allowing  us  to  write  the  general  electric  field  in  region  1  as  Ey  =  AelklZ  +  Be~lkiz  and  in 
region  3  as  Ey  =  Celk?,L  +  De~lksL  where  ki  =  ^  with  c o  being  the  frequency  and  c  the 
speed  of  light.  To  obtain  the  general  solution  in  region  2  we  must  solve  the  sourceless, 
macroscopic  Maxwell’s  equations  with  a  spatially  varying  dielectric  function.  This  produces 
three  uncoupled  scalar  wave  equations,  one  for  each  component  of  the  electric  field  in  the 
graded  region.  However,  the  fact  that  E  =  yEy  in  regions  1  and  3  along  with  the  boundary 
conditions  and  the  divergence  condition  V  •  E  =  —  A—  (A  Ez  will  force  Ex  =  Ez  =  0  inside 
region  2.  We  are  left  with  the  scalar  wave  equation  for  Ey  given  by 


d2Ey  e(z)^ 
dz 2  +  c2  y 


0. 


(160) 
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The  general  solution  for  Ey  when  the  index  of  refraction  varies  according  to  \J e  (z)  = 
n  (z)  —  ni  +  mz  is 


EV  =  C1 


n(z)  T  f  Tm2  (z) 

u 


+  C<2 


fnn2(z) 


m  4  \  A  m  J  *  V  m  4  ^  A  m 

where  Ji  and  Yi  are  quarter-order  Bessel  functions  of  the  hrst  and  second  kind,  respec- 

4  4 

tively  (see  later  details  for  geometric  series  expressions  for  Ji  and  Yi).  Introducing  the 
shorthand  notation 


(161) 
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we  can  write  the  characteristic  matrix  M  (ref. 
the  held  En  and  its  ^-derivative  from  z  =  0  to  z  = 


nn2  (z)  A 
A  m  ) 

TTU2  (z)  A 

A  m  J 
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(163) 


Born  and  Wolf  [58]),  which  propagates 
L,  as 


M  = _ 1 _  (  h  (L)  cf  (0)  -  h!  (0)  g  (L)  h  (0)  g  (L)  -  h  (L)  g  (0)  \ 

h  (0)  g'  (0)  -  h'  (0)  9  (0)  l  h!  (L)  g'  (0)  -  h!  (0)  g'  (L)  h  (0)  g'  ( L )  -h'(L)g  (0)  J 

(164) 

where  primes  denote  differentiation  with  respect  to 

Because  the  index  is  continuous  at  the  grade  boundaries,  the  E  and  B  fields  must  be 
continuous  there.  For  our  solution  this  requires  the  continuity  of  Ey  and  its  derivative. 
Applying  the  boundary  conditions  and  using  M  to  propagate  Ey  and  its  derivative  across 
the  grade,  we  can  relate  the  fields  in  regions  1  and  3  according  to 


/  C  ik3L  +  De-ik3L  \  /  A  +  B  \ 

=  Ml  .  (165) 

\  ik3  (Ceik3L  -  De~ik3L)  J  [ik^A-B)  J 

First  we  solve  for  transmission  and  reflectance  coefficients  for  top-down  incidence  where 
A  =  1  and  D  =  0.  This  gives 


Cteik3L  _  —2ki 

IgifciO  _mn£3  _  +  i  ( mi2kik3  -  m2 1)  ’ 

Bte~lk i°  _  mnA:3  -  m22ki  +  i  (m^iAq  +  m21) 
IgifciO  — fflnfc3  -  m22ki  +  i  ( mnkiks  -  m2 1) 


(166) 

(167) 
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where  the  rriij  are  the  elements  of  M  and  we  have  used  the  fact  that  M  is  unimodular. 
Now  we  solve  for  the  coefficients  with  bottom-up  incidence  where  D  —  1  and  A  =  0 
obtaining 


Bbe-ikl°  _  2k3 

le-iksL  muk3  +  m22ki  -  i  [mi2kik3  -  m2 1)  ’ 
ChelkaL  _  mnk3  -  m22ki  -  i  (mi2kik3  +  m2i) 

1  e~iksL  muk3  +  m22ki  -  i  ( mi2kik3  -  m2 1) ' 


(168) 

(169) 


Note  that  k i  ty,  =  k3tt  which  implies  that  in  the  absence  of  absorption  the  transmitted 
phase  is  independent  of  propagation  direction.  From  these  relations,  we  obtain  the  desired 
phase  angles  a  and  (3  by  expressing  the  coefficients  in  exponential  form  as: 


tt  =  tmeiai,  rt  =  rtoel131,  tb  =  tboeia 2,  rt  =  rboe102 .  (170) 

Figure  57  plots  the  various  phase  angles  as  a  function  of  grade  length  for  several  popular 
material  systems. 

Once  the  phase  angles  for  the  desired  grade  thickness  are  known,  the  thicknesses  of  the 
constant  index  layers  can  be  chosen  using  the  phase  matching  conditions  (equations  (154) 
and  (156)).  The  reflectivity  spectrum  of  the  complete  structure  can  be  calculated  using  the 
transfer  matrix  technique  by  using  equation  (164)  for  the  graded  layer. 

For  completeness,  the  fractional  Bessel  functions  of  the  first  and  second  kind  can  be 
evaluated  using  the  following  expressons: 


k\  p( 

k= o  v 


(— l)fc(f)*'+2fe 

V  +  k  +  1) 


,  Yv(x)  = 


Ju(x)  COs(7 TU)  —  J-U{x ) 
sin(7Ti') 


(171) 


The  geometric  series  for  Ji 


rn2(z) 

Am 


converges  within  calculable  error  after  ten  or  twenty 


terms  for  reasonable  grade  lengths  and  index  values. 


GROWTH,  FABRICATION,  AND  CHARACTERIZATION  OF  LOW-THRESHOLD 
VCSELS 

In  this  section  is  detailed  both  the  enhancements  in  core  capability  as  well  as  results  of 
our  work  in  the  creation  of  very-low  threshold  VCSELs.  As  mentioned  in  the  introduction, 
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FIG.  57:  Phase  angle  plots  for  several  popular  materials,  a)  AlAs/GaAs  grade  at  980  nm,  b) 
AlAs/Alo.i6Gao.84As  at  850  nm,  c)  AlAs/Alo.3Gao.7As  at  780  nm,  d)  AlSb/Alo.iGao.gSb  at  1550 
nm 

much  of  the  initial  effort  of  this  project  went  into  two  main  areas-improvements  to  modelling 
and  simulations  for  correct  designs,  and  then  improvement  to  growth  capabilities,  which  is 
the  heart  of  any  successful  device-focused  research. 

Growth  Capability  Development 

Our  initial  forays  into  developing  ultra-low  threshold  VCSELs  began  with  growth  of  cal¬ 
ibration  structures,  such  as  multiple  quantum  well  (MQW)  samples  for  determination  of 
emission  wavelength  and  photoluminescence  (PL)  intensity,  and  n—  and  p-doped  DBR  mir¬ 
rors.  The  mirror  growths  served  not  only  to  calibrate  the  growth  rates,  when  accompanied 
by  fits  of  the  accompanying  reflection  spectra,  but  also  as  a  means  of  studying  the  conduction 
of  current  through  the  respective  doped  mirror  structure. 
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Initially,  our  aim  was  to  “walk  before  you  run,”  and  we  took  the  approach  of  attempting 
to  demonstrate  lasing  in  a  conventional  etched  post,  pump-through-the- mirrors  VCSEL  de¬ 
sign  before  proceeding  into  state-of-the-art,  very-low  threshold,  oxide-apertured,  intracavity 
contacted  designs. 

During  these  initial  stages,  many  growth  issues  and  challenges  facing  the  current  MBE 
configuration  became  clear.  For  instance,  the  growth  ovens  and  controllers  were  relatively 
unstable  over  the  long  times  required  for  a  VCSEL  growth  (>  14  hours  in  typical  cases). 
Indeed,  reflection  high-energy  electron  diffraction  (RHEED)  measurements  taken  prior  to 
and  after  VCSEL  growth  showed  significant  discrepancies.  At  first,  we  attempted  to  correct 
this  problem  by  adding  optical  monitors  to  the  MBE  system,  and  using  fits  to  reflectance 
data  at  various  stages  in  sample  growth  to  correct  for  such  growth  rate  drifts.  The  Varian 
Gen  II  MBE  system  used  in  this  project  is  not  readily  adaptable  for  in-situ  real-time 
monitoring  of  growth.  However,  this  system  does  have  affixed  to  it  a  sample  preparation 
chamber  used  to  load  multiple  wafers  onto  a  wheeled  substrate  holder  prior  to  an  interlock 
connecting  to  the  actual  growth  chamber.  This  intermediate  chamber  has  optical  viewports 
allowing  the  potential  for  reflection  (R( A))  and  photoluminescence  (PL( A))  measurements 
while  under  high- vacuum  conditions.  In  this  way,  we  sought  to  break  the  growth  of  a 
VCSEL  into  stages-bottom  mirror  growth,  growth  of  active  region,  and  top  mirror  growtli  - 
and  perform  such  quasi-in-situ  PL( A)  and  R( A)  measurements.  The  information  from  these 
spectra  and  resultant  fits  could  then  allow  growth  rate  re-calibration  without  the  need  to 
load  new  RHEED  samples  into  the  chamber,  and  also  without  the  need  to  break  high- 
vacuum.  Although  this  is  not  as  useful  as  true  real-time  growth  rate  monitoring,  it  does 
significantly  shorten  the  turnaround  time  and  offer  the  opportunity  for  increased  success 
for  each  growth  run.  Figure  58  shows  an  example  of  R( A)  and  PL( A)  spectra  measured 
for  such  an  application.  This  figures  shows  that  the  fit  values  for  Al„Gai_xAs  and  GaAs 
growth  rates  are  98.2%  and  97.8%  of  the  ideal  rates,  respectively.  Also,  this  figure  shows 
how  the  photoluminescence  spectrum  peaks  near  971  nm,  whereas  the  cavity  dip  measured 
after  growth  of  the  bottom  mirror  and  active  region  of  the  VCSEL  indicates  a  reflectivity 
“notch”  at  967.3  nm,  indicating  relatively  good  alignment  between  gain  peak  and  cavity 
resonance. 

Although  this  ability  of  quasi-in-situ  growth  monitoring  helped  immensely,  it  was  still  far 
from  allowing  us  to  generate  a  working  VCSEL.  Two  other  significant  enhancements  to  the 
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FIG.  58:  Measured  reflectivity  and  photoluminescence  spectra  for  a  VCSEL  sample  at  various 
stages  during  growth.  The  red  curve  shows  the  calculated  fit  to  only  the  bottom  DBR  spectrum 
(not  shown);  the  blue  curve  is  R( A)  for  the  DBR  +  active  region,  with  an  indicated  cavity  resonance 
at  967.3  nm;  the  magenta  curve  shows  PL( A)  at  this  same  stage  in  growth,  with  a  peak  at  971.0 
nm,  indicating  suitable  alignment  of  the  gain  to  the  cavity. 

MBE  tool  were  then  performed.  The  first  was  a  complete  overhaul  of  the  control  code  and 
machinery  interfacing  to  the  machine.  In  a  sweeping  re- vamp  of  this  system,  we  replaced  not 
only  a  host  of  convention  growth  effusion  cells  with  larger  capacity  (400  g)  SUMO  cells  from 
EPI  MBE  Products.  The  larger  capacity  cells  not  only  would  allow  longer  operation  of  the 
machine  between  re-filling  of  the  cells,  but  were  also  demonstrated  to  be  more  stable  for  long¬ 
term  growth  operations.  Along  with  the  cells  came  a  host  of  upgrades  for  the  temperature 
controllers,  and  a  new  interfacing  system  (a  National  Instruments  FieldPoint  Module )  to 
allow  faster  communications  between  a  desktop  computer  system  to  all  of  the  various  controls 
of  the  MBE  (including  oven  temperatures  settings,  shutter  actuators,  and  substrate  rotation 
control).  The  underlying  software  architecture  for  device  growth  was  then  separated  into 
two  elements  of  code.  The  first  part  was  a  National  Instruments  Lab  View  interface  control 
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starttable 

vcsel.def 

delay(100,as) 

repeat  41  times 

696  galas(si)  “bottom  mirror” 

812.69  alga2as(si)  “bottom  mirror” 

end 

608.74  gal  as(si)  “oxide  phase  match  layer” 

180  digitalgrade(galas(si),alga2as(si),steps=9)  “grade  up” 

superlattice(10,alga2as(si),40.18,alas(si),periods=13)  “98%  layer” 

180  di gitalgrade(alga2as(si), gal  as(si), step  s=9)  “grade  down” 

318.74  galas  “note:  setup  carbon” 

set(si,850) 

delay(220,as)  “cool  down,  arsenic  up” 

set(subst,510) 

100  galas 

repeat  3  times 

delay(17,in,gal,as)  “QW” 

100  galas  “barrier” 

end 

delay(120,as)  “arsenic  down” 

set(subst,600) 

318.74  galas 

180  digitalgrade(galas(c),alga2as(c),steps=9)  “grade  up” 

superlattice(10,alga2as(c),40.1 8,alas(c),periods=13)  “98%  layer” 

180  di gitalgrade(alga2as(c), gal  as(c), step s=9)  “grade  down” 

500  galas(c) 

repeat  31  times 

180  digitalgrade  (gal  as(c),alga2as(c),steps=9) 

672.1  alga2as(c) 

180  digitalgrade  (alga2as(c),galas(c),steps=9) 

520.57  galas(c) 

end 

180  digitalgrade  (gal  as(c),alga2as(c),steps=9) 

672.1  alga2as(c) 

180  digitalgrade  (alga2as(c),galas(c),steps=9) 

408.3  galas(c) 

200  galas(c) 

delay(600,as) 

set(subst,100) 

endtableQ 

FIG.  59:  Sample  Epitaxial  Description  Language  input  deck  for  the  instructions  for  growth  of  a 
980-nm  VCSEL  via  MBE. 

code  to  perform  all  of  the  necessary  hardware  operations,  including  “graceful”  shutdown  of 
the  system  for  maintenance.  The  second  piece  of  control  software,  the  so-called  Epitaxial 
Description  Language  (EDL),  was  a  Microsoft  Visual  C++  compiler  code  that  renders  text 
input  from  a  Microsoft  Word  table  into  a  series  of  commands  capable  of  being  read  by  the 
Lab  View  code  and  summarily  translated  into  machine  operable  instructions.  This  overhaul 
of  the  control  software  for  the  MBE  system,  although  time  consuming  and  laborious,  resulted 
in  an  elegant  and  flexible  “language”  for  growth  of  various  structures  by  this  system.  An 
example  of  an  EDL  compiler  code  table  for  the  growth  of  an  entire  VCSEL  is  demonstrated 
in  Fig.  59.  Examination  of  this  input  deck  will  show  the  inclusion  of  digital  alloy  grading  in 
both  the  top  p-doped  DBR  mirror,  as  well  as  to  form  high  Al  mole  fraction  oxide  aperture 
layers. 

The  second  major  modification  to  this  system  was  the  inclusion  of  a  carbon  tetrabro- 
mide  (CBrzt)  p-doping  source.  It  has  been  found  that  the  use  of  carbon  as  a  p-doping 
source  provides  two  main  benefits:  (1)  the  ability  to  dope  at  much  higher  concentrations 
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FIG.  60:  Plot  of  activated  dopant  concentration  versus  depth  in  sample  for  a  carbon-doped  GaAs 
wafer. 

(~  1020cm-3),  and  (2)  the  fact  that  carbon  as  a  p-type  dopant  does  not  diffuse  nearly  as 
much  as  more  typical  dopants,  such  as  beryllium.  In  all,  these  combine  to  allow  better  tai¬ 
lored  p-dopant  profiles.  As  this  flavor  of  dopant  generally  presents  more  material  optical  loss 
than  n-dopants,  one  can  readily  see  the  advantages  of  being  able  to  place  heavily  p-doped 
layers  in  regions  of  the  optical  cavity  where  the  lasing  mode  field  intensity  is  highest,  and 
vice  versa.  Figure  60  shows  a  plot  of  activated  dopant  concentration  versus  depth  in  sample 
for  a  p-GaAs  epitaxial  growth.  This  was  performed  as  a  calibration  sample  to  verify  both 
our  ultimate  doping  levels  with  this  new  oven,  and  also  the  control  with  valve  position  that 
this  source  enables. 

Processing  Capability  Development 

The  key  processing  steps  for  most  semiconductor  devices,  after  successful  epitaxial 
growth,  are  typically  (1)  Deposition  of  metal  contacts,  (2)  ability  to  etch  various  mate¬ 
rial  with  high  fidelity,  and  (3)  ability  to  isolate  or  passivate  materials  for  profiled  conduction 
of  current  through  a  device.  This  is  equally  true  for  VCSEL  processing,  but  in  our  case,  an¬ 
other  key  processing  step  includes  the  ability  to  perform  selective  oxidation  in  high  aluminum 
mole  fraction  AffGa^j-As  materials,  the  benefits  of  which  has  been  previously  discussed. 
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FIG.  61:  SEM  micrograph  of  an  ICP/RIE  etch  of  an  AlAs-GaAs  etalon  structure.  The  individual 
layers  comprising  the  top  (8  periods)  and  bottom  (8.5  periods)  DBRs  are  clearly  visible.  Of  note 
is  the  sidewall  smoothness  and  anisotropic  nature  of  the  etch. 

Ohmic  contacts  to  standard  III-V  materials  are  well  established,  and  we  have  made 
no  effort  to  improve  upon  the  voluminous  amount  of  prior  existing  work.  The  reader  may 
therefore  reference  any  standard  text  on  III-V  material  processing  for  suggested  metal  recipes 
for  VCSEL  fabrication.  For  our  effort,  p-metal  deposition  was  typically  Ti/Au  or  Ti/Pt/Au, 
whereas  n-metal  depositions  were  of  the  Ni/Ge/Au/Ni/Au  variety. 

For  mesa  etching,  however,  we  employed  a  (then)  state-of-the-art  technology,  inductively- 
coupled-plasm  reactive  ion  etching  (ICP/RIE),  to  perform  mesa  isolation  and  VCSEL  device 
definition.  The  tool  used  was  a  Plasma-Therm  (now  Unaxis )  model  SLR770  ICP/RIE 
system.  This  system  is  design  etching  a  single  wafer  at  a  time,  and  is  a  loadlock  system  with 
6  process  gases.  Our  etch  recipe  used  30  seem  of  BCR  with  10  seem  of  CR,  and  the  etch  was 
done  typically  under  the  following  tool  conditions:  pressure  of  4  mTorr,  wafer  temperature 
of  10  0  C,  utilizing  600  W  of  ICP  power  and  50  W  of  RIE  power.  Results  from  etching 
an  “empty  etalon”  structure  consisting  of  8  periods  of  alternating  top  AlAs/GaAs  DBR 
layers  and  8.5  periods  of  bottom  DBR  layers  with  a  1-A  GaAs  etalon  cavity  layer  are  shown 
in  Fig.  61.  From  this  SEM  micrograph,  it  is  evident  that  we  are  obtaining  very  smooth 
and  vertical  sidewalls  using  this  etch  recipe.  For  sake  of  comparison,  each  alternating  light 
(AlAs)  and  dark  (GaAs)  layer  is  comprising  either  mirror  stack  is  nominally  789  A  or  696 
A  thick,  respectively. 

Furthermore,  in  order  to  enable  us  to  precisely  control  the  etch  depth,  we  outfitted 
this  system  with  a  homemade  reflectance  etch  monitor  system  (Fig.  62).  This  system  is 
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FIG.  62:  Schematic  cross  section  of  the  in-house  ICP/RIE  system,  complete  with  custom-built 
reflectance  monitoring  system. 

comprised  of  a  laser  diode  (nominally  A iaser  ~  789  nm ),  a  silicon  photodiode  detector,  an 
imaging  camera  with  laser  line  filter  (to  block  all  but  a  small  portion  of  the  laser  diode 
signal,  to  prevent  saturation  during  alignment),  and  a  removable  beamsplitter.  The  output 
of  the  photodiode  is  then  sent  to  a  desktop  PC  to  collect  the  reflected  signal  as  a  function 
of  time.  The  etch  trace  collected  by  this  system  for  the  sample  shown  in  Fig.  61  is  given  in 
Fig.  63. 

As  an  early  demonstration  of  the  powerful  capability  of  this  system,  we  processed  VC- 
SEL  material  grown  by  the  Army  Research  Laboratory,  Sensors  and  Electron  Devices  Direc¬ 
torate  (at  the  time,  under  the  direction  of  George  Simonis)  .  This  ARL  structure,  denoted 
ARL1236 ,  was  an  intracavity  contacted  device  incorporating  dual  oxide  apertures  on  either 
side  of  the  gain  region.  Figure  64  shows  our  results  in  etching  this  sample.  In  this  figure, 
the  a;- axis  represents  depth  into  the  sample  in  /a in.  with  the  vertical  axis  plotting  either  the 
refractive  index  profile  of  the  device,  or  the  etch  trace  obtained  from  our  etching  system. 
The  shaded  areas  represent  the  regions  of  top-  or  bottom-contact  areas  for  this  device,  nom¬ 
inally  |A  thick,  with  A  the  design  wavelength  (980  nm)  and  n  the  refractive  index  of  the 
material.  Also  shown  on  this  plot  is  the  final  etch  profile  of  this  device  taken  on  a  Tencor 
P-10  stylus  profilometer.  This  figure  demonstrates  our  ability  to  etch  device  with  the  tight 
tolerances  necessary  to  form  intracavity  contact  VCSELs. 

With  our  ability  to  etch  samples  and  deposit  Ohmic  contacts,  the  final  piece  of  the 
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FIG.  63:  Plot  of  reflected  signal  versus  time  for  the  etched  sample  of  Figure  61. 
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FIG.  64:  Plot  of  reflected  signal  versus  time  for  an  intracavity  contacted  VCSEL  grown  by  the 
Army  Research  Laboratory. 
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FIG.  65:  Schematic  diagram  of  our  custom-built  oxidation  furnace  complete  with  capability  for 
in-situ  monitoring. 

puzzle  for  fabrication  of  oxide-apertured  VCSELs  was  the  formation  of  the  oxide  apertures 
themselves.  In  order  to  reproducibly  define  small  features  by  oxidizing  AlGaAs  layers,  it 
is  essential  to  have  good  control  over  the  oxidation  reaction.  To  push  VCSEL  structures 
to  their  ultimate  threshold  limits,  these  apertures  must  constrict  to  1  pm  and  below.  [59]. 
Unfortunately,  the  oxidation  reaction  in  a  conventional  furnace  is  hard  to  control,  making  it 
difficult  to  reproducibly  define  such  small  features.  We  have  integrated  a  glass  viewport  into 
a  cold-walled  oxidation  chamber  to  enable  in-situ  optical  monitoring  of  the  sample  during 
oxidation  (see  Fig.  65).  To  prevent  water  vapor  from  condensing  on  the  cold  window  and 
obscuring  the  sample,  the  reaction  chamber  is  operated  at  low  pressure  (5  Torr).  This  low 
pressure  also  allows  rapid  evacuation  of  the  chamber,  which  immediately  halts  the  reaction. 
To  gain  additional  control,  we  reduced  the  oxidation  temperature  to  325  °C,  consequently 
slowing  the  oxidation  rate  to  2  pm  per  hour.  The  combination  of  optical  monitoring,  slow 
oxidation  rate,  and  rapid  shutoff  greatly  increases  our  control  over  the  reaction  and  allows 
us  to  define  small  features  precisely. 


90 


Real-time  in-situ  optical  measurements  of  AlAs  oxidation  rates  were  performed  using  this 
system  and  the  results  were  compared  with  a  standard  model  [60].  Oxide-semiconductor 
distributed  Bragg  reflectors  were  also  fabricated  and  measured,  yielding  highly-reflective 
mirrors  suitable  for  vertical-cavity  surface-emitting  laser  fabrication. 

To  verify  that  these  oxides  have  suitable  optical  and  mechanical  properties  for  DBR  mirror 
fabrication,  we  also  measured  the  reflectance  of  oxidized  DBR  features.  The  low  oxidation 
rate  makes  it  impractical  to  oxidize  large  features.  Therefore  we  constructed  an  imaging 
system  capable  of  measuring  the  reflectance  of  small  (~20  //in)  features.  The  samples 
under  study  were  composed  of  five  periods  of  GaAs  (695  A)  and  AlAs  (1530  A)  grown 
by  molecular  beam  epitaxy  on  a  GaAs  substrate.  Vertical  posts,  with  square  and  circular 
cross  sections  ranging  in  size  from  5  to  40  microns,  were  reactive-ion-etched  into  the  wafer 
using  the  following  procedure.  The  wafer  was  first  covered  by  SiaN.j  using  plasma-enhanced 
chemical  vapor  deposition  and  then  coated  with  photoresist.  The  resist  was  patterned  using 
conventional  lithographic  techniques,  the  SiaNi  was  anisotropically  etched  in  a  CF3H/O2 
plasma,  and  the  remaining  resist  was  removed.  The  patterned  SisN4  was  used  as  a  mask  for 
an  anisotropic  BCly/He  plasma  etch  of  the  epitaxial  layers,  and  was  subsequently  removed 
by  isotropic  etching  in  a  CF4/O2  plasma.  Finally,  the  exposed  AlAs  layers  in  the  wafer 
were  protected  from  atmospheric  water  vapor  by  re-depositing  500  Aof  SRN4  over  the  entire 
wafer.  The  resulting  features  exhibited  side-wall  angles  less  than  5°  from  perpendicular. 
All  oxidation  samples  were  prepared  from  this  single  wafer:  immediately  prior  to  oxidation, 
samples  were  cleaved  from  the  wafer  and  the  SI3N4  protective  layer  was  removed  by  a  CF4/O2 
plasma  etch. 

The  oxidation  furnace,  schematically  depicted  in  Figure  65,  comprises  a  conventional 
sample  heater  with  an  integral  thermocouple  to  control  the  sample  temperature.  Pure 
water  vapor  is  delivered  to  the  chamber  at  a  fixed  flow  rate  through  a  vapor-source  mass- 
flow  controller.  A  feedback  control  system  maintains  the  chamber  pressure  at  5  Torr:  the 
system  consists  of  a  capacitance-manometer  pressure  sensor  driving  a  butterfly  throttle- 
valve  on  a  mechanical  vacuum  pump.  As  the  AlAs  layers  are  oxidized,  the  near-infrared 
reflectivity  of  the  DBR  samples  changes  dramatically.  By  illuminating  the  sample  with 
near-infrared  light  and  viewing  the  reflection  with  a  silicon  CCD  camera  attached  to  a  long- 
working-distance  microscope,  we  easily  discriminate  between  the  oxidized  and  unoxidized 
regions.  The  dramatic  contrast  is  evident  in  the  near-infrared  photograph  shown  in  Figure 
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FIG.  66:  Near  infrared  image  of  a  sample  during  the  oxidation  process.  The  high  refractive  index 
contrast  between  the  unoxidized  AlAs  (n  ~  3.0),  and  native  oxide  (n  ~  1.6)  is  readily  apparent  in 
this  image. 

66. 

Samples  were  oxidized  in  the  system  at  temperatures  ranging  from  325  to  400  °C.  The 
color  of  the  oxidized  samples  was  identical  for  all  oxidation  conditions.  Thick  AlAs  samples 
that  were  oxidized  downward  from  the  top  surface  of  the  wafer  were  found  to  be  somewhat 
unstable  mechanically,  tending  to  delaminate  from  the  substrate  material  just  at  the  point 
when  the  layer  was  fully  oxidized.  However,  laterally  oxidized  samples  showed  no  tendencies 
to  delaminate  under  any  of  the  oxidation  conditions  discussed  here,  and  in  fact  were  found 
to  be  quite  rugged.  In  addition,  unprotected  AlAs  samples  were  found  to  degrade  quite 
quickly  (in  a  matter  of  a  few  days)  in  ambient  conditions,  but  oxidized  samples,  even  those 
oxidized  at  low  temperatures,  show  no  signs  of  degradation  after  nearly  one  year.  Since  the 
purpose  of  this  study  is  to  develop  a  process  to  fabricate  VCSELs  using  laterally  oxidized 
layers,  we  hereafter  consider  only  laterally  oxidized  samples. 

By  viewing  the  samples  through  the  CCD  camera  and  charting  the  lateral  extent  of  the 
oxide  as  a  function  of  time,  we  measured  the  oxidation  rate  without  removing  the  sample 
from  the  chamber  or  interrupting  the  reaction.  The  results  are  plotted  in  Figure  67.  The 
temperature  range  studied  was  limited  at  the  low  end  by  the  length  of  time  needed  to 
reliably  measure  an  appreciable  oxidation  distance,  and  was  limited  at  the  high  end  by  the 
components  of  the  oxidation  system,  notably  the  elastomer  seals  on  the  chamber.  For  some 
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FIG.  67:  Oxidation  distance  as  a  function  of  time  for  several  temperatures.  All  oxidations  were 
carried  out  at  a  water  vapor  pressure  of  5  Torr  and  a  flow  rate  of  500  seem. 

samples,  a  small  time  delay  (approximately  5-15  min)  elapsed  before  noticeable  oxidation 
occurred;  the  delay  did  not  vary  systematically  with  oxidation  temperature.  This  delay 
could  result  from  uncontrolled  surface  oxidation  in  the  brief  time  between  SRN4  removal 
and  sample  loading. 

VCSEL  Device  Results 

As  a  result  of  all  of  the  previous  underlying  work,  we  designed,  grew,  and  fabricated  a  las¬ 
ing  structure  aimed  at  sub-mA  threshold  current  operation.  The  resultant  structure  design 
was  a  dual  oxide  apertured  structure  employing  21  carbon-doped  p-type  GaAs/Alo.gGao.iAs 
DBR  mirror  layers  (complete  with  180  A  digital  alloy  grades  between  the  materials),  a  1-A 
cavity  region  with  3  80-A  In0.2Ga0.sAs  as  the  active  region.  The  bottom  DBR  consisted  of 
41  pairs  of  silicon  doped  GaAs/Alo.gGao.iAs  n-type  mirror  layers  without  any  grading  be¬ 
tween  the  layers.  As  seen  in  Figures  68,  69,  we  have  made  an  extra-cavity  contacted  device 
(for  ease  of  fabrication)  that,  when  operating  at  room  temperature,  displays  a  threshold 
current  of  ith  =  122  /iA  with  a  threshold  voltage  of  Vth  =  2.76  V.  As  detailed  in  Figure 
68,  we  see  the  following  behavior  for  lasing  wavelength,  A/ase,  as  a  function  of  the  oxide 
aperture  radius,  pox'-  \ase  holds  steady  for  devices  with  larger  oxide  apertures  (little  effect 
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FIG.  68:  Plot  of  ith  and  A iase  as  a  function  of  oxide  aperture  radius  for  an  in-house  design,  grown, 
and  fabricated  YCSEL  structure.  See  the  accompanying  text  for  the  device  description. 

on  the  microcavity  resonance  condition);  as  the  size  of  the  aperture  shrinks,  the  wavelength 
actually  redshifts  slightly  (presumably  due  to  device  heating  due  to  increased  current  crowd¬ 
ing  effects);  and  eventually  the  lasing  emission  blueshifts,  as  predicted  by  both  the  VFEM 
and  WIMP  simulations.  As  our  “cold  cavity”  simulations  were  not  performed  under  a  cou¬ 
pled  optical  +  thermal  +  electrical  phenomenology,  the  effects  of  current  crowding  inducing 
heating  in  these  structures  was  not  incorporated. 

One  of  the  critical  features  for  obtaining  very  efficient,  low-threshold  VCSELs  is  to  ensure 
that  the  gain  peak  aligns  with  the  cavity  resonance  under  anticipated  operating  conditions. 
Indeed,  this  was  part  of  the  underlying  motivation  for  development  of  the  WIMP  and  VFEM 
tools-to  incorporate  such  shifts  a  priori  in  our  designs.  One  is  never  guaranteed,  however, 
that  growth  will  match  the  design  conditions.  To  test  that  the  gain  peak  was  close  to 
the  cavity  resonance  for  the  smallest  oxide  aperture  devices  (first  data  points,  lowest  ith  in 
Figures  68,  69),  we  performed  similar  experiments  in  a  temperature  controlled  probe-stand. 
It  is  well  known  that  the  gain  peak  shifts  much  more  rapidly  with  temperature  than  the 
cavity  resonance  (see,  e.g.,  [61]  and  references  therein).  Hence,  but  performing  measurements 
of  ith  as  a  function  of  temperature,  we  may  infer  the  amount  of  mismatch  in  gain  and  cavity 
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FIG.  69:  Plot  of  jth  and  Vth  as  a  function  of  oxide  aperture  radius  for  an  in-house  design,  grown, 
and  fabricated  VCSEL  structure.  See  the  accompanying  text  for  the  device  description. 
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FIG.  70:  Plot  of  V^,  and  A iase  as  a  function  of  temperature  for  one  of  the  oxide  aperture 
VCSELs  shown  in  Figure  68.  The  minimum  for  ith  and  Vth  occurs  near  17  °C,  indicating  good 
match  between  design  and  the  ultimate  device  results. 

resonances.  Figure  70  shows  the  results  of  this  experiment. 
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From  this  figure,  it  is  evident  that  our  design  was  very  close  to  ideal,  in  than,  for  a 
structure  with  pox  =  1-8/ an,  the  minima  in  both  ith  and  Vfi,  occurs  just  slightly  below  room 
temperature,  nominally  17  °C.  For  completeness,  the  shift  in  A iase  with  temperature  is  also 
plotted. 

SUMMARY  AND  FUTURE  EFFORTS 

In  the  AFOSR  supported  in-house  research  project,  we  embarked  on  an  effort  to  generate 
extremely  low-threshold  VCSELs.  To  support  this  effort,  state  of  the  art  simulation  suites 
were  developed  to  outline  the  optical  physics  of  microcavity  structures,  including  etched 
post  designs,  and  oxide  apertures  layers  with  or  without  spatial  tapered  profiles.  The  two 
methods,  the  weighted  index  method  with  parasitic  mode  coupling  (WIMP),  and  the  full 
vector  finite-element  method  (VFEM),  has  been  applied  to  many  VCSEL  structures.  We 
have  succeeded  in  predicting  lasing  wavelength  blueshifts  with  decreasing  oxide  aperture 
diameters,  both  for  published  results  and  for  our  own  in-house  made  devices.  We  have 
successfully  managed  to  describe  the  relative  design  merits  for  inclusion  of  oxide  apertures 
in  microcavity  resonator  designs,  again  through  explaining  and  interpreting  the  heretofore 
confusing  results  found  in  the  literature.  We  have  also  shown  the  benefits  toward  applying 
tapered  oxide  aperture  designs  for  “blunting”  the  diffraction  losses  incurred  by  “blunt”  aper¬ 
tures.  Finally,  in  the  design  arena,  we  have  developed  quasi- analytic  means  for  calculating 
the  transmission  and  reflection  properties  of  graded-layer  structures,  and  showed  the  utility 
of  this  generalized  interface  approach  for  design  of  DBR  mirrors  and  VCSELs  using  digital 
alloying  grading  at  the  material  interfaces. 

Our  in-house  growth  and  fabrication  yielded  VCSELs  operating  cw  at  room  temperature 
with  threshold  currents  as  low  as  122  //A.  For  a  manufacturable  (i.e.  high-yield)  structure 
utilizing  current  injection  through  the  p-  and  n-DBR  mirrors,  this  was  less  than  a  factor  of 
2  from  the  best  published  results  at  the  time,  and  less  than  a  factor  of  14  from  the  best 
devices  period. 

As  a  result  of  this  research  effort,  we  have  also  greatly  increase  the  growth  and  fabrication 
capability  at  AFRL.  An  overhaul  of  the  Varian  GEN  II  MBE  machine  used  for  this  project 
has  results  in  VCSEL  class  material.  We  have  also  demonstrated  a  custom  built  steam 
oxidation  furnace  for  creation  of  native  oxides  of  aluminum  from  AlGaAs  material,  com- 
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FIG.  71:  Plot  the  reduced  voltage  needed  to  drive  a  bipolar  cascade  edge  emitting  laser.  A  detailed 
description  of  this  device  is  not  in  the  scope  of  this  report,  and  the  reader  is  encourage  to  consult 
reference  [62]  for  further  details  on  this  newer  effort. 

plete  with  in-situ  monitoring  allowing  for  high-precision  control  of  oxide-apertured  VCSEL 
structures. 

Future  efforts  for  optoelectronics  and  semiconductor  device  research  will  see  spillover 
benefits  from  the  efforts  made  on  this  project.  Indeed,  we  have  embarked  on  efforts  to 
generate  interband  (so-called  “bipolar”)  cascade  lasers  for  the  purpose  of  radiofrequency- 
photonic  links.  These  devices,  which  incorporate  Esaki  tunnel  junctions  as  small-resistance 
interconnects  between  lasing  active  regions,  are  predicted  to  have  superior  noise  and  overall 
link  gain  characteristics  compared  to  tradition  laser  diodes,  a  feature  we  aim  to  demonstrate 
experimentally.  Indeed,  we  have  already  shown  single-  and  multi-color  emission  from  edge 
emitting  structures  (see  Figure  71  and  [62],  for  example),  and  we  are  now  pursuing  VCSEL 
geometries  for  device  improvements.  This,  and  other  efforts,  will  be  the  subject  of  a  future 
report. 
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